*COMDECK GRIDA
      COMMON/GRIDA/X(170,85),Y(170,85)
*COMDECK GRIDB
      COMMON/GRIDB/DS(170),THETA(170),AAA(170),BBB(170),DSOB(170),
     1  THETOB(170),CCC(170),DDD(170),P1(170),Q1(170),R1(170),S1(170),
     2  OBANGS(170),DIST(241),AIRFX(241),AIRFY(241),CAMBRX(170),
     3  CAMBRY(170),XOB(170),YOB(170),XMIN(100),XMAX(100),YMIN(100),
     4  YMAX(100),JLE1A(2),JLE2A(2),NORDA(2),MAXITA(2),JCORN(4)
*COMDECK GRIDC
      COMMON/GRIDC/MAXIT,JS,JE,JINC,JTEBOT,JTETOP,KS,KE,KINC,JMAX,KMAX,
     1  NTETYP,JPRT,NLETYP,JN1,JN2,NPLT,NIBDST,NOBSHP,NOBDST,NSEQ,
     2  JCAMBR,NAIRF,NORD,JAIRF,JDIST,NOUT,NDS,OMEGA,OMEGP,OMEGQ,
     3  OMEGR,OMEGS,TODEG,PLIM,QLIM,RLIM,SLIM,RADOB,XOBCNT,XLEFT,
     4  XRIGHT,YBOTOM,YTOP,RCORN,PI,TWOPI,PIO2,ALAMF,ALAMR,XLE,XTE,
     5  TEOPEN,BINN,TR,XTFRAC,ROTANG,ROTCTR,DSI,AAAI,BBBI,CCCI,DDDI,
     6  DSOBI,THETAI,THOBI,WAKEP,NOBCAS
*COMDECK GRID1
      NAMELIST/GRID1/JMAX,KMAX,NTETYP,NAIRF,JAIRF,NIBDST,NDS,DSI,JTEBOT,
     1  JTETOP,TR,XLE,XTE,NOBSHP,NOBDST,RADOB,XLEFT,XRIGHT,YBOTOM,YTOP,
     2  RCORN,ALAMF,ALAMR,NORDA,MAXITA,JPRT,NPLT,NOUT
*COMDECK GRID2
      NAMELIST/GRID2/NLETYP,BINN,JDIST,JCAMBR,XTFRAC,ROTANG,ROTCTR,
     1  XOBCNT,OMEGA,OMEGP,OMEGQ,OMEGR,OMEGS,PLIM,QLIM,RLIM,SLIM,
     2  DSOBI,THETAI,THOBI,AAAI,BBBI,CCCI,DDDI,TEOPEN,WAKEP,NOBCAS
*COMDECK GRID3
      NAMELIST/GRID3/DS,DSOB,THETA,THETOB,DIST,AIRFX,AIRFY,CAMBRX,
     1  CAMBRY,OBANGS,XOB,YOB,XMIN,XMAX,YMIN,YMAX,AAA,BBB,CCC,DDD
*DECK BLOCK
      OVERLAY(GRAPE,0,0)
C
      BLOCK DATA
C
*CALL GRIDA
C
*CALL GRIDB
C
*CALL GRIDC
C
      DATA JMAX/  100/,KMAX/  49/,NTETYP/  1/,NAIRF/  2/,JAIRF/  0/,
     1  NIBDST/  1/,JTEBOT/  15/,JTETOP/  86/,TR/  .12/,XLE/  0./,
     2  XTE/  1./,NOBSHP/  1/,NOBDST/  1/,RADOB/  6./,XLEFT/  -6./,
     3  XRIGHT/  6./,  YBOTOM/  -4./,YTOP/  4./,RCORN/  1./,ALAMF/  0./,
     4  ALAMR/  0./,NORDA/  4,  1/,MAXITA/  200,  100/,JPRT/  0/,
     5  NPLT/  0/,NOUT/  1/,NOBCAS/  0/
C
      DATA NLETYP/  1/,BINN/  1.1/,JDIST/  0/,JCAMBR/  0/,XTFRAC/  1./,
     1  ROTANG/  0./,ROTCTR/  0./,XOBCNT/  0./,OMEGA/  1.3/,
     2  OMEGP,OMEGQ,OMEGR,OMEGS/  4*.3/,  PLIM,QLIM,RLIM,SLIM/  4*1./,
     3  NDS/  2/,DSI/  0./,DSOBI/  0./,THETAI/  0./,THOBI/  0./,
     4  AAAI,BBBI,CCCI,DDDI/  4*0./,TEOPEN/  0./,WAKEP/  1./
C
      DATA DS/  170*.01/,DSOB/  170*.2/,THETA/  170*90./,
     1  THETOB/  170*90./,DIST/  241*0./,AIRFX/  241*0./,
     2  AIRFY/  241*0./,CAMBRX/  170*0./,CAMBRY/  170*0./,
     3  OBANGS/  170*0./,XOB/  170*0./,YOB/  170*0./,XMIN,XMAX,YMIN,
     4  YMAX/  100*0.,  100*0.,  100*0.,  100*0./,AAA,BBB,CCC,DDD/
     5  170*.45,  170*.45,  170*.45,  170*.45/
C
      DATA JCORN/4*0/
C
      END
*DECK MAIN00
      PROGRAM MAIN00(INPUT,TAPE5=INPUT,OUTPUT,TAPE6=OUTPUT,TAPE7,TAPE15)
C
*CALL GRIDA
C
*CALL GRIDB
C
*CALL GRIDC
C
C        GENERATE AND STORE THE GRID
      CALL OVERLAY(5HGRAPE,1,0)
C
C        PLOT THE GRID
      IF(NPLT.GT.0)  CALL OVERLAY(5HGRAPE,2,0)
C
      STOP
C
      END
*DECK MAIN10
      OVERLAY(1,0)
C
      PROGRAM MAIN10
C
      COMMON/SCRACH/DUM(241,6)
C
      CALL INPUT
C
      CALL INCHK
C
      CALL INNER
C
      CALL OUTER
C
      CALL SOLVE
C
      CALL OUTPUT
C
      RETURN
C
      END
*DECK CKSMTH
      SUBROUTINE CKSMTH(F,TOL,JMAX,J,KEY)
C
C        THIS SUBROUTINE CHECKS FOR SMOOTHNESS IN ARRAY F BY PASSING
C        A PARABOLA THROUGH THE NEAREST 3 POINTS AND OBSERVING
C        THE DIFFERENCE BETWEEN THE ACTUAL VALUE FOR F(J) AND THE
C        CORRESPONDING POINT ON THE PARABOLA.
C
      DIMENSION F(4)
C
      IF(J.GT.1)  GO TO 1
      F1=F(2)
      F2=F(3)
      F3=F(4)
      X1=2.
      X2=3.
      X3=4.
      XX=1.
      GO TO 2
C
    1 IF(J.GT.2)  GO TO 3
      F1=F(1)
      F2=F(3)
      F3=F(4)
      X1=1.
      X2=3.
      X3=4.
      XX=2.
      GO TO 2
C
    3 IF(J.LT.JMAX)  GO TO 4
      F1=F(JMAX-3)
      F2=F(JMAX-2)
      F3=F(JMAX-1)
      X1=1.
      X2=2.
      X3=3.
      XX=4.
C
    2 X1SQ=X1**2
      X2SQ=X2**2
      X3SQ=X3**2
      D1=X2SQ*X3-X3SQ*X2
      D2=X1SQ*X3-X3SQ*X1
      D3=X1SQ*X2-X2SQ*X1
      RD=1./(D1-D2+D3)
      A=X3*F2-X2*F3-X3*F1+X1*F3+X2*F1-X1*F2
      B=X2SQ*F3-X3SQ*F2-X1SQ*F3+X3SQ*F1+X1SQ*F2-X2SQ*F1
      C=F1*D1-F2*D2+F3*D3
      GO TO 5
C
    4 F1=F(J-2)
      F2=F(J-1)
      F3=F(J+1)
      XX=3.
C
      RD=-1./6.
      A=-2.*F1+3.*F2-F3
      B=12.*F1-15.*F2+3.*F3
      C=-16.*F1+12.*F2-2.*F3
C
    5 FF=(A*XX**2+B*XX+C)*RD
      DIF=TOL*ABS(FF)
      KEY=0
      IF(DIF.GT.0..AND.ABS(FF-F(J)).GT.DIF)  KEY=1
      RETURN
      END
*DECK CSPLIN
      SUBROUTINE CSPLIN(XX,YY,YYP,X,Y,N1,N2,J1,J2,KYYP)
      DIMENSION XX(1),YY(1),YYP(1),X(1),Y(1)
      COMMON/SCRACH/A(241),B(241),C(241),D(241),F(241),H(241)
C
C          CUBIC SPLINE INTERPOLATION
C  X,Y ARRAYS ARE TO BE INTERPOLATED
C  YY ARE FOUND INTERPOLATES CORRESPONDING TO XX
C  J1,,J2, ARE INDICE LIMITS ON X,Y
C  N1, N2 ARE INDICE LIMITS ON XX ( ALSO YY)
C  DIMENSION OF ARRAYS CARRIED IN FROM OUTSIDE, X(J) MUST BE MONOTONE
C FORMULA FROM NUMERICAL METHODS BY DAHLQUIST,BJORCK,ANDERSON
C                                                 JLS FEB. 77
C        MODIFIED JUNE, 79 BY RLS.  THIS SUBROUTINE WAS SET UP TO
C        PRINT AN ERROR MESSAGE AND QUIT IF ANY ELEMENT OF XX WAS
C        OUTSIDE OF THE RANGE GIVEN BY X.  THE PROBLEM WAS THAT ROUND-
C        OFF COULD CAUSE THIS CONDITION TO INADVERTANTLY BE VIOLATED.
C        THE MODIFICATION WAS TO RELAX THE CONDITION SLIGHTLY, I.E.,
C        TO ALLOW A LITTLE BIT OF EXTRAPOLATION.
C
C        MODIFIED AGAIN OCTOBER, 79 BY RLS TO CALCULATE SLOPE OF YY.
C
C      FIRST FIND DERIVATIVE LIKE TERMS THAT ARE COEFFICIENTS
      JA = J1 + 1
      JB = J2-1
      DO 10 J=JA,J2
      H(J) = X(J) - X(J-1)
   10 D(J) =(Y(J) - Y(J-1))/H(J)
      DO 12 J=JA,JB
      A(J) = H(J+1)
      B(J) = 2.*( H(J) + H(J+1))
   12 F(J) = 3.*( H(J)*D(J+1) + H(J+1)*D(J))
      B(J1) = 2.
      H(J1) = 1.
      F(J1) = 3.*D(JA)
      A(J2) = 1.
      B(J2) = 2.
      F(J2) = 3.*D(J2)
      CALL TRIB(A,B,H,C,F,J1,J2,1)
C
C  INTERPOLATION , X(J) ARRAY MUST BE MONOTONE
      FUDGE=1.E-10*ABS(X(J2)-X(J1))
      J = J1
      I = J1 + 1
      DO 20 N=N1,N2
   28 IF(X(J)-FUDGE.LE.XX(N).AND.X(I)+FUDGE.GE.XX(N))  GO TO 30
      IF( X(I) - XX(N))22,22,24
   22 J = J +1
      I = I+1
      IF( I .GT. J2) GO TO 27
      GO TO 28
   24 J = J -1
      I = I-1
      IF( J.LT. J1) GO TO 27
      GO TO 28
   27 WRITE(6,600)
  600 FORMAT(55H0FATAL ERROR IN CSPLIN.  VALUE AT WHICH WE ARE SUPPOSED,
     1  52H TO INTERPOLATE IS OUTSIDE OF GIVEN TABLE OF VALUES.)
      STOP
   30 T = (XX(N) - X(J))/H(I)
      TT = 1. - T
      YY(N) = T*Y(I) + TT*Y(J) + H(I)*T*TT*( (F(J) - D(I))*TT
     1   -(F(I)-D(I))*T)
      IF(KYYP.EQ.1)  YYP(N)=F(J)+2.0*(-2.0*F(J)-F(I)+3.0*D(I))*T
     1  +3.*(F(J)+F(I)-2.0*D(I))*T**2
   20 CONTINUE
      RETURN
      END
*DECK IC
      SUBROUTINE IC
C
C        THIS SUBROUTINE FILLS THE X, Y, P1, Q1, R1, AND S1 ARRAYS
C        WITH INITIAL CONDITIONS.  X AND Y ARE FILLED USING A
C        PRE-DETERMINED EXPONENTIAL STRECHING.
C
      COMMON/SCRACH/DISTIC(85),DUM(1361)
C
*CALL GRIDA
C
*CALL GRIDB
C
*CALL GRIDC
C
      DATA EPS,DSM,NCALL/.2,.01,0/
C
      NCALL=NCALL+1
      K2=KS+KINC
      K3=KE-KINC
      JF=JS
      JL=JE
      IF(NTETYP.LT.3)  GO TO 4
      JF=1
      JL=JMAX
C
    4 DISTIC(KS)=0.
      DO 1 K=K2,KE,KINC
      KK=K-KINC
    1 DISTIC(K)=DISTIC(KK)+DSM*(1.+EPS)**(K-KS-KINC)
      RLAST=1./DISTIC(KE)
      IF(NCALL.EQ.1.AND.JPRT.GE.0)  WRITE(6,100)
      DO 2 K=KS,KE,KINC
      DISTIC(K)=DISTIC(K)*RLAST
      IF(NCALL.EQ.1.AND.JPRT.GE.0)  WRITE(6,101)  K,DISTIC(K)
    2 CONTINUE
C
      DO 3 J=JF,JL,JINC
      DX=X(J,KE)-X(J,KS)
      DY=Y(J,KE)-Y(J,KS)
      P1(J)=0.
      Q1(J)=0.
      R1(J)=0.
      S1(J)=0.
C
      DO 3 K=K2,K3,KINC
      X(J,K)=X(J,KS)+DISTIC(K)*DX
    3 Y(J,K)=Y(J,KS)+DISTIC(K)*DY
      RETURN
  100 FORMAT(54H0INITIAL CONDITIONS WILL BE ACCORDING TO THE FOLLOWING,
     1  23H EXPONENTIAL STRECHING:/3X,1HK,9X,6HDISTIC)
  101 FORMAT(1X,I3,5X,F10.6)
      END
*DECK INCHK
      SUBROUTINE INCHK
C
*CALL GRIDB
C
*CALL GRIDC
C
      LOGICAL FATAL
C
*CALL GRID1
C
*CALL GRID2
C
*CALL GRID3
C
      DATA FATAL/.FALSE./
      DATA TOLDS/  .05/,TOLDSO/  .05/, TOLTH/  .1/,TOLTHO/  .1/,
     1  TOLDST/  .05/,TOLAFX/  .1/,TOLAFY/  .1/,TOLCBX/  .5/,
     2  TOLCBY/  .5/,TOLOBA/  .1/,TOLAAA,TOLBBB,TOLCCC,TOLDDD/  4*.05/
      DATA TOLXOB/  .1/,TOLYOB/  .1/
C
C        SET CONSTANTS INVOLVING PI.
      PI=3.1415926535
      TWOPI=2.*PI
      PIO2=PI*.5
      TODEG=180./PI
      TORAD=PI/180.
C
C        PRINT BANNER PAGE AND INPUT VARIABLES.
      IF(JPRT.LT.0)  GO TO 275
      WRITE(6,376)
      WRITE(6,377)
      WRITE(6,GRID1)
      WRITE(6,GRID2)
  275 CONTINUE
C
C        *************************************************************
C        *                                                           *
C        *     CHECK VARIABLES IN NAMELIST /GRID1/.                  *
C        *                                                           *
C        *************************************************************
C
      IF(JMAX.GE.4.AND.JMAX.LE.170)  GO TO 2
      WRITE(6,101)
      FATAL=.TRUE.
C
    2 IF(KMAX.GE.4.AND.KMAX.LE.85)  GO TO 3
      WRITE(6,102)
      FATAL=.TRUE.
C
    3 IF(NTETYP.GE.1.AND.NTETYP.LE.3)  GO TO 4
      WRITE(6,103)
      FATAL=.TRUE.
C
    4 IF(NAIRF.GE.1.AND.NAIRF.LE.5)  GO TO 1
      WRITE(6,104)
      FATAL=.TRUE.
C
    1 IF(JAIRF.GE.0.AND.JAIRF.LE.241)  GO TO 5
      WRITE(6,100)
      FATAL=.TRUE.
C
    5 IF(NIBDST.GE.1.AND.NIBDST.LE.5)  GO TO 6
      WRITE(6,105)
      FATAL=.TRUE.
C
    6 IF(NIBDST.NE.3.AND.NIBDST.NE.5.OR.NAIRF.EQ.5)  GO TO 26
      WRITE(6,126)
      FATAL=.TRUE.
C
   26 IF(NAIRF.LT.5.OR.JAIRF.GE.4)  GO TO 28
      WRITE(6,125)
      FATAL=.TRUE.
C
   28 IF(TR.GE.0..AND.TR.LT.1.)  GO TO 7
      WRITE(6,106)
      FATAL=.TRUE.
C
    7 CONTINUE
C
    8 IF(NIBDST.EQ.5)  GO TO 9
      IF(XTE.GT.XLE)  GO TO 9
      WRITE(6,108)
      FATAL=.TRUE.
C
    9 IF(NOBSHP.GE.1.AND.NOBSHP.LE.6)  GO TO 40
      WRITE(6,109)
      FATAL=.TRUE.
   40 IF(NTETYP.LE.2.AND.NOBSHP.LE.3.OR.
     1   NTETYP.LE.2.AND.NOBSHP.EQ.6.OR.
     2   NTETYP.EQ.3.AND.NOBSHP.GE.4)  GO TO 42
      WRITE(6,140)
      FATAL=.TRUE.
   42 IF(NOBSHP.LE.2.OR.NOBSHP.EQ.4.OR.NOBDST.EQ.1)  GO TO 10
      WRITE(6,142)
      FATAL=.TRUE.
C
   10 IF(NOBDST.GE.1.AND.NOBDST.LE.3)  GO TO 11
      WRITE(6,110)
      FATAL=.TRUE.
C
   11 CONTINUE
C
   12 IF(NIBDST.EQ.5.OR.NOBSHP.EQ.1)  GO TO 13
      IF(XLEFT.LT.XLE)  GO TO 13
      WRITE(6,112)
      FATAL=.TRUE.
C
   13 IF(NIBDST.EQ.5.OR.NOBSHP.EQ.1)  GO TO 14
      IF(XRIGHT.GT.XLE)  GO TO 14
      WRITE(6,113)
      FATAL=.TRUE.
C
   14 IF(YBOTOM.LT.0.)  GO TO 15
      WRITE(6,114)
      FATAL=.TRUE.
C
   15 IF(YTOP.GT.0.)  GO TO 16
      WRITE(6,115)
      FATAL=.TRUE.
C
   16 IF(RCORN.GT.0..AND.RCORN.LE..5*(YTOP-YBOTOM).AND.RCORN.LE..5*
     1  (XRIGHT-XLEFT))  GO TO 17
      WRITE(6,116)
      FATAL=.TRUE.
C
   17 IF(ABS(ALAMF).LT.90.)  GO TO 18
      WRITE(6,117)
      FATAL=.TRUE.
C
   18 IF(ABS(ALAMR).LT.90.)  GO TO 43
      WRITE(6,118)
      FATAL=.TRUE.
C
   43 IF(NOBSHP.EQ.3.OR.NOBSHP.EQ.5)  GO TO 19
      IF(ALAMR.EQ.0.)  GO TO 37
      WRITE(6,143)
      FATAL=.TRUE.
      GO TO 37
   19 IF(ALAMF.NE.0..OR.ALAMR.NE.0.)  GO TO 37
      WRITE(6,137)
      FATAL=.TRUE.
C
   37 IF(NORDA(1).GE.0)  GO TO 20
      WRITE(6,119)
      FATAL=.TRUE.
C
   20 IF(NORDA(2).GE.0)  GO TO 21
      WRITE(6,120)
      FATAL=.TRUE.
C
   21 IF(MAXITA(1).GE.0)  GO TO 22
      WRITE(6,121)
      FATAL=.TRUE.
C
   22 IF(MAXITA(2).GE.0.)  GO TO 23
      WRITE(6,122)
      FATAL=.TRUE.
C
   23 CONTINUE
C
   24 IF(NPLT.GE.0.AND.NPLT.LE.100)  GO TO 25
      WRITE(6,124)
      FATAL=.TRUE.
C
   25 IF(NOUT.GE.0.AND.NOUT.LE.3)  GO TO 93
      WRITE(6,193)
      FATAL=.TRUE.
C
C        SET OR CHECK JTEBOT AND JTETOP.
   93 IF(NTETYP.EQ.3)  GO TO 27
      JTEBOT=1
      JTETOP=JMAX
      IF(JPRT.GE.0)  WRITE(6,127)
      GO TO 29
   27 IF(JTEBOT.GE.1.AND.JTETOP.GT.JTEBOT.AND.JTETOP.LE.160.AND.
     1  JTEBOT-1.EQ.JMAX-JTETOP)  GO TO 29
      WRITE(6,128)
      FATAL=.TRUE.
C
C        DECIDE WHETHER OR NOT SEQUENCING CAN BE USED.
   29 IF(MAXITA(1).GT.0)  GO TO 94
      WRITE(6,194)
      GO TO 35
   94 IF(((KMAX-1)/3)*3.EQ.KMAX-1)  GO TO 95
      WRITE(6,195)
      GO TO 35
   95 GO TO (31,32,33),NTETYP
   31 IF((JMAX/4)*4.EQ.JMAX)  GO TO 34
      WRITE(6,129)
      GO TO 35
   32 IF((JMAX/3)*3.EQ.JMAX)  GO TO 34
      WRITE(6,130)
      GO TO 35
   33 IF(((JMAX-1)/3)*3.EQ.JMAX-1)  GO TO 34
      WRITE(6,131)
   35 NSEQ=1
      MAXITA(2)=MAXITA(2)+MAXITA(1)
      NORDA(2)=NORDA(2)+NORDA(1)
      WRITE(6,378)  MAXITA(2),NORDA(2)
      GO TO 36
   34 IF(JPRT.GE.0)  WRITE(6,132)
      NSEQ=0
   36 CONTINUE
C
C        *************************************************************
C        *                                                           *
C        *     CHECK VARIABLES IN NAMELIST /GRID2/.                  *
C        *                                                           *
C        *************************************************************
C
      IF(NLETYP.GE.1.AND.NLETYP.LE.3)  GO TO 52
      WRITE(6,146)
      FATAL=.TRUE.
C
   52 IF(NLETYP.EQ.1.OR.NIBDST.EQ.5)  GO TO 254
      WRITE(6,353)
      FATAL=.TRUE.
C
  254 IF(NLETYP.EQ.1.OR.JMAX.EQ.JAIRF)  GO TO 280
      WRITE(6,380)
      FATAL=.TRUE.
C
  280 IF(BINN.GT.1)  GO TO 53
      WRITE(6,152)
      FATAL=.TRUE.
C
   53 IF(NIBDST.NE.4.OR.JDIST.GE.4.AND.JDIST.LE.241)  GO TO 30
      WRITE(6,133)
      FATAL=.TRUE.
C
   30 IF(JCAMBR.EQ.0.OR.JCAMBR.GE.4.AND.JCAMBR.LE.170)  GO TO 54
      WRITE(6,153)
      FATAL=.TRUE.
C
   54 IF(XTFRAC.GT.0.)  GO TO 55
      WRITE(6,154)
      FATAL=.TRUE.
C
   55 IF(ABS(ROTANG).LT.90.)  GO TO 56
      WRITE(6,155)
      FATAL=.TRUE.
C
   56 CONTINUE
C
   57 IF(NOBSHP.GT.1.OR.XOBCNT+RADOB.GT.XTE.AND.XOBCNT-RADOB.LT.XLE)
     1  GO TO 58
      WRITE(6,111)
      FATAL=.TRUE.
C
   58 IF(OMEGA.GT.0..AND.OMEGA.LT.2.)  GO TO 59
      WRITE(6,158)
      FATAL=.TRUE.
C
   59 IF(OMEGP.GE.0..AND.OMEGP.LT.2.)  GO TO 60
      WRITE(6,159)
      FATAL=.TRUE.
C
   60 IF(OMEGQ.GE.0..AND.OMEGQ.LT.2.)  GO TO 61
      WRITE(6,160)
      FATAL=.TRUE.
C
   61 IF(OMEGR.GE.0..AND.OMEGR.LT.2.)  GO TO 62
      WRITE(6,161)
      FATAL=.TRUE.
C
   62 IF(OMEGS.GE.0..AND.OMEGS.LT.2.)  GO TO 63
      WRITE(6,162)
      FATAL=.TRUE.
C
   63 IF(PLIM.GT.0..AND.PLIM.LT.100.)  GO TO 64
      WRITE(6,163)
      FATAL=.TRUE.
C
   64 IF(QLIM.GT.0..AND.QLIM.LT.100.)  GO TO 65
      WRITE(6,164)
      FATAL=.TRUE.
C
   65 IF(RLIM.GT.0..AND.RLIM.LT.100.)  GO TO 66
      WRITE(6,165)
      FATAL=.TRUE.
C
   66 IF(SLIM.GT.0..AND.SLIM.LT.100.)  GO TO 67
      WRITE(6,166)
      FATAL=.TRUE.
C
   67 IF(NDS.EQ.1.OR.NDS.EQ.2)  GO TO 68
      WRITE(6,167)
      FATAL=.TRUE.
C
   68 IF(DSI.GE.0.)  GO TO 69
      WRITE(6,168)
      FATAL=.TRUE.
C
   69 IF(DSOBI.GE.0.)  GO TO 70
      WRITE(6,169)
      FATAL=.TRUE.
C
   70 IF(ABS(THETAI-90.).LE.90..OR.NIBDST.EQ.5)  GO TO 71
      WRITE(6,170)
      FATAL=.TRUE.
C
   71 IF(ABS(THOBI-90.).LE.90..OR.NOBSHP.EQ.6)  GO TO 72
      WRITE(6,171)
      FATAL=.TRUE.
C
   72 IF(AAAI.GE.0.)  GO TO 73
      WRITE(6,172)
      FATAL=.TRUE.
C
   73 IF(BBBI.GE.0.)  GO TO 74
      WRITE(6,173)
      FATAL=.TRUE.
C
   74 IF(CCCI.GE.0.)  GO TO 75
      WRITE(6,174)
      FATAL=.TRUE.
C
   75 IF(DDDI.GE.0.)  GO TO 76
      WRITE(6,175)
      FATAL=.TRUE.
C
   76 IF(TEOPEN.GE.0.)  GO TO 38
      WRITE(6,176)
      FATAL=.TRUE.
C
   38 IF(WAKEP.GT.0..AND.WAKEP.LT.2.)  GO TO 39
      WRITE(6,139)
      FATAL=.TRUE.
   39 CONTINUE
      IF(NOBCAS.EQ.0.OR.NOBCAS.EQ.1)  GO TO 41
      WRITE(6,141)
      FATAL=.TRUE.
   41 CONTINUE
C
C        FILL SOME ARRAYS.
      IF(NDS.EQ.1.OR.DSI.EQ.0.)  GO TO 77
      DO 78 J=1,JMAX
   78 DS(J)=DSI
C
   77 IF(DSOBI.EQ.0.)  GO TO 79
      DO 80 J=1,JMAX
   80 DSOB(J)=DSOBI
C
   79 IF(THETAI.EQ.0.)  GO TO 81
      DO 82 J=1,JMAX
   82 THETA(J)=THETAI
C
   81 IF(THOBI.EQ.0.)  GO TO 83
      DO 84 J=1,JMAX
   84 THETOB(J) = THOBI
C
   83 IF(AAAI.EQ.0.)  GO TO 85
      DO 86 J=1,JMAX
   86 AAA(J)=AAAI
C
   85 IF(BBBI.EQ.0.)  GO TO 87
      DO 88 J=1,JMAX
   88 BBB(J)=BBBI
C
   87 IF(CCCI.EQ.0.)  GO TO 89
      DO 90 J=1,JMAX
   90 CCC(J)=CCCI
C
   89 IF(DDDI.EQ.0.)  GO TO 91
      DO 92 J=1,JMAX
   92 DDD(J)=DDDI
   91 CONTINUE
C
C        *************************************************************
C        *                                                           *
C        *     CHECK VARIABLES IN NAMELIST /GRID3/.                  *
C        *                                                           *
C        *************************************************************
C
C        WRITE INPUT VALUES IN GRID3
      IF(JPRT.GE.0.)  WRITE(6,GRID3)
C
C        CHECK ARRAY DS FOR RANGE AND SMOOTHNESS.
      NWRT=0
      DO 221 J=1,JMAX
      IF(DS(J).GT.0.)  GO TO 221
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,321)  J,DS(J)
      FATAL=.TRUE.
  221 CONTINUE
C
      NWRT=0
      DO 201 J=1,JMAX
      CALL CKSMTH(DS    ,TOLDS,JMAX,J,KEY)
      IF(KEY.EQ.0)  GO TO 201
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,301)  J,DS    (J)
  201 CONTINUE
C
C        CHECK ARRAY DSOB FOR RANGE AND SMOOTHNESS.
  200 NWRT=0
      DO 222 J=1,JMAX
      IF(DSOB(J).GT.0.)  GO TO 222
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,322)  J,DSOB(J)
      FATAL=.TRUE.
  222 CONTINUE
C
      NWRT=0
      DO 202 J=1,JMAX
      CALL CKSMTH(DSOB  ,TOLDSO,JMAX,J,KEY)
      IF(KEY.EQ.0)  GO TO 202
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,302)  J,DSOB  (J)
  202 CONTINUE
C
C        CHECK ARRAY THETA FOR RANGE AND SMOOTHNESS.
      IF(NIBDST.EQ.5)  GO TO 410
      NWRT=0
      DO 223 J=1,JMAX
      IF(ABS(THETA(J)-90.).LT.90.)  GO TO 223
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,323)  J,THETA(J)
      FATAL=.TRUE.
  223 CONTINUE
C
  410 NWRT=0
      DO 203 J=1,JMAX
      CALL CKSMTH(THETA ,TOLTH,JMAX,J,KEY)
      IF(KEY.EQ.0)  GO TO 203
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,303)  J,THETA (J)
  203 CONTINUE
C
C        CHECK ARRAY THETOB FOR RANGE AND SMOOTHNESS.
      IF(NOBSHP.EQ.6)  GO TO 411
      NWRT=0
      DO 224 J=1,JMAX
      IF(ABS(THETOB(J)-90.).LT.90.)  GO TO 224
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,324)  J,THETOB(J)
      FATAL=.TRUE.
  224 CONTINUE
C
  411 NWRT=0
      DO 204 J=1,JMAX
      CALL CKSMTH(THETOB,TOLTHO,JMAX,J,KEY)
      IF(KEY.EQ.0)  GO TO 204
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,304)  J,THETOB(J)
  204 CONTINUE
C
C        CHECK ARRAY DIST FOR RANGE, SMOOTHNESS, AND MONOTONICITY.
      IF(NIBDST.NE.4)  GO TO 239
      NWRT=0
      DO 225 J=1,JDIST
      IF(J.EQ.1.AND.DIST(J).EQ.0.)  GO TO 225
      IF(J.EQ.JDIST.AND.DIST(J).EQ.1.)  GO TO 225
      IF(J.GT.1.AND.J.LT.JDIST.AND.DIST(J).GT.0..AND.DIST(J).LT.1.)
     1  GO TO 225
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,325)  J,DIST  (J)
      FATAL=.TRUE.
  225 CONTINUE
C
      NWRT=0
      DO 205 J=1,JDIST
      CALL CKSMTH(DIST,TOLDST,JDIST,J,KEY)
      IF(KEY.EQ.0)  GO TO 205
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,305)  J,DIST  (J)
  205 CONTINUE
C
      NWRT=0
      DO 240 J=2,JDIST
      JM=J-1
      IF(DIST(J).GT.DIST(JM))  GO TO 240
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,340)  J,JM
      FATAL=.TRUE.
  240 CONTINUE
C
C        CHECK ARRAY AIRFX FOR SMOOTHNESS.  ALSO MAKE SURE THAT THE
C        ELEMENTS OF AIRFX ARE NOT ALL THE SAME.
  239 IF(JAIRF.LT.4)  GO TO 215
      NWRT=0
      DO 206 J=1,JAIRF
      CALL CKSMTH(AIRFX,TOLAFX,JAIRF,J,KEY)
      IF(KEY.EQ.0)  GO TO 206
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,306)  J,AIRFX (J)
  206 CONTINUE
C
      DO 244 J=2,JAIRF
      IF(AIRFX(J).NE.AIRFX(J-1))  GO TO 215
  244 CONTINUE
      WRITE(6,344)
      FATAL=.TRUE.
C
C        IF AIRFOIL IS TO BE INTERPOLATED, CHECK THAT IT GOES ALL THE
C        WAY FROM TRAILING-EDGE AROUND TO TRAILING EDGE.
      IF(NIBDST.EQ.5)  GO TO 215
      IF(AIRFX(1).EQ.AIRFX(JAIRF))  GO TO 215
      WRITE(6,379)
      FATAL=.TRUE.
C
C        CHECK ARRAY AIRFY FOR SMOOTHNESS.
  215 IF(JAIRF.LT.4)  GO TO 216
      NWRT=0
      DO 207 J=1,JAIRF
      CALL CKSMTH(AIRFY,TOLAFY,JAIRF,J,KEY)
      IF(KEY.EQ.0)  GO TO 207
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,307)  J,AIRFY (J)
  207 CONTINUE
C
C        CHECK ARRAY CAMBRX FOR RANGE, SMOOTHNESS, AND MONOTONICITY.
  216 IF(JCAMBR.LT.4)  GO TO 217
      NWRT=0
      DO 228 J=1,JCAMBR
      IF(J.EQ.1.AND.CAMBRX(J).EQ.0.)  GO TO 228
      IF(J.EQ.JCAMBR.AND.CAMBRX(J).EQ.1.)  GO TO 228
      IF(J.GT.1.AND.J.LT.JCAMBR.AND.CAMBRX(J).GT.0..AND.CAMBRX(J).
     1  LT.1.)  GO TO 228
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,328)  J,CAMBRX(J)
      FATAL=.TRUE.
  228 CONTINUE
C
      NWRT=0
      DO 208 J=1,JCAMBR
      CALL CKSMTH(CAMBRX,TOLCBX,JMAX,J,KEY)
      IF(KEY.EQ.0)  GO TO 208
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,308)  J,CAMBRX(J)
  208 CONTINUE
C
      NWRT=0
      DO 241 J=2,JCAMBR
      JM=J-1
      IF(CAMBRX(J).GT.CAMBRX(JM))  GO TO 241
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,341)  J,JM
      FATAL=.TRUE.
  241 CONTINUE
C
C        CHECK ARRAY CAMBRY FOR SMOOTHNESS.
  217 IF(JCAMBR.LT.4)  GO TO 218
      NWRT=0
      DO 209 J=1,JCAMBR
      CALL CKSMTH(CAMBRY,TOLCBY,JMAX,J,KEY)
      IF(KEY.EQ.0)  GO TO 209
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,309)  J,CAMBRY(J)
  209 CONTINUE
C
C        CHECK ARRAY OBANGS FOR RANGE, SMOOTHNESS, AND MONOTONICITY.
  218 IF(NOBDST.NE.3)  GO TO 235
      IF(NOBSHP.NE.4)  GO TO 236
      OBAL=-ATAN2(YBOTOM,XRIGHT)*TODEG
      OBAH=360.-ATAN2(YTOP,XRIGHT)*TODEG
      OBAL=FLOAT(IFIX(OBAL*100.+.5))*.01
      OBAH=FLOAT(IFIX(OBAH*100.+.5))*.01
      NWRT=0
      DO 230 J=1,JMAX
      IF(J.EQ.1.AND.ABS(OBANGS(J)-OBAL).LT..00001)  GO TO 230
      IF(J.EQ.JMAX.AND.ABS(OBANGS(J)-OBAH).LT..00001)  GO TO 230
      IF(J.GT.1.AND.J.LT.JMAX.AND.OBANGS(J).GT.OBAL.AND.OBANGS(J).LT.
     1  OBAH)  GO TO 230
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,330)  J,OBANGS(J)
      FATAL=.TRUE.
  230 CONTINUE
      IF(NWRT.EQ.1)  WRITE(6,343)  OBAL,OBAH
C
  236 NWRT=0
      DO 210 J=1,JMAX
      CALL CKSMTH(OBANGS,TOLOBA,JMAX,J,KEY)
      IF(KEY.EQ.0)  GO TO 210
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,310)  J,OBANGS(J)
  210 CONTINUE
C
      NWRT=0
      DO 242 J=2,JMAX
      JM=J-1
      IF(OBANGS(J).GT.OBANGS(JM))  GO TO 242
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,342)  J,JM
      FATAL=.TRUE.
  242 CONTINUE
C
C        CHECK ARRAY XOB FOR SMOOTHNESS
  235 IF(NOBSHP.NE.6)  GO TO 270
      NWRT=0
      DO 273 J=1,JMAX
      CALL CKSMTH(XOB,TOLXOB,JMAX,J,KEY)
      IF(KEY.EQ.0)  GO TO 273
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,373)  J,XOB(J)
  273 CONTINUE
C
C        CHECK ARRAY YOB FOR SMOOTHNESS
  270 IF(NOBSHP.NE.6)  GO TO 272
      NWRT=0
      DO 274 J=1,JMAX
      CALL CKSMTH(YOB,TOLYOB,JMAX,J,KEY)
      IF(KEY.EQ.0)  GO TO 274
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,374)  J,YOB(J)
  274 CONTINUE
C
C        CHECK ARRAY AAA FOR RANGE AND SMOOTHNESS.
  272 NWRT=0
      DO 231 J=1,JMAX
      IF(AAA(J).GE.0.)  GO TO 231
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,331)  J,AAA   (J)
      FATAL=.TRUE.
  231 CONTINUE
C
      NWRT=0
      DO 211 J=1,JMAX
      CALL CKSMTH(AAA   ,TOLAAA,JMAX,J,KEY)
      IF(KEY.EQ.0)  GO TO 211
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,311)  J,AAA   (J)
  211 CONTINUE
C
C        CHECK ARRAY BBB FOR RANGE AND SMOOTHNESS.
      NWRT=0
      DO 232 J=1,JMAX
      IF(BBB(J).GE.0.)  GO TO 232
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,332)  J,BBB   (J)
      FATAL=.TRUE.
  232 CONTINUE
C
      NWRT=0
      DO 212 J=1,JMAX
      CALL CKSMTH(BBB   ,TOLBBB,JMAX,J,KEY)
      IF(KEY.EQ.0)  GO TO 212
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,312)  J,BBB   (J)
  212 CONTINUE
C
C        CHECK ARRAY CCC FOR RANGE AND SMOOTHNESS.
      NWRT=0
      DO 233 J=1,JMAX
      IF(CCC(J).GE.0.)  GO TO 233
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,333)  J,CCC   (J)
      FATAL=.TRUE.
  233 CONTINUE
C
      NWRT=0
      DO 213 J=1,JMAX
      CALL CKSMTH(CCC   ,TOLCCC,JMAX,J,KEY)
      IF(KEY.EQ.0)  GO TO 213
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,313)  J,CCC   (J)
  213 CONTINUE
C
C        CHECK ARRAY DDD FOR RANGE AND SMOOTHNESS.
      NWRT=0
      DO 234 J=1,JMAX
      IF(DDD(J).GE.0.)  GO TO 234
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,334)  J,DDD   (J)
      FATAL=.TRUE.
  234 CONTINUE
C
      NWRT=0
      DO 214 J=1,JMAX
      CALL CKSMTH(DDD   ,TOLDDD,JMAX,J,KEY)
      IF(KEY.EQ.0)  GO TO 214
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,314)  J,DDD   (J)
  214 CONTINUE
C
C        CHECK ARRAYS XMIN, XMAX, YMIN, AND YMAX.
      IF(NPLT.EQ.0)  GO TO 219
      NWRT=0
      DO 238 N=1,NPLT
      IF(XMAX(N).GT.XMIN(N).OR.YMAX(N).GT.YMIN(N))  GO TO 238
      IF(NWRT.EQ.0)  WRITE(6,300)
      NWRT=1
      WRITE(6,316)  XMIN(N),XMAX(N),YMIN(N),YMAX(N),N
  238 CONTINUE
  219 CONTINUE
C
C        IN THE CASE OF NIBDST=5, THE USER HAS SUPPLIED AN AIRFOIL, AND
C        THOSE AIRFOIL POINTS ARE TO BE USED EXACTLY AS READ.  IN THAT
C        CASE, REDEFINE XLE AND XTE TO CORRESPOND TO THE GIVEN AIRFOIL.
C        THEN DO ALL OF THE TESTS USING XLE AND XTE.
      IF(NIBDST.NE.5)  GO TO 252
      XLE=AIRFX(1)
      XTE=AIRFX(1)
      DO 44 J=2,JAIRF
      IF(AIRFX(J).LT.XLE)  XLE=AIRFX(J)
      IF(AIRFX(J).GT.XTE)  XTE=AIRFX(J)
   44 CONTINUE
      WRITE(6,144)  XLE,XTE
      IF(NOBSHP.EQ.1.OR.XLEFT.LT.XLE)  GO TO 250
      WRITE(6,112)
      FATAL=.TRUE.
  250 IF(NOBSHP.EQ.1.OR.NOBSHP.EQ.6.OR.XRIGHT.GT.XTE)  GO TO 251
      WRITE(6,113)
      FATAL=.TRUE.
  251 IF(NOBSHP.GT.1.OR.XOBCNT+RADOB.GT.XTE.AND.XOBCNT-RADOB.LT.XLE)
     1  GO TO 252
      WRITE(6,111)
      FATAL=.TRUE.
C
C        CONVERT ALL OF THE ANGLES TO RADIANS.
  252 ALAMF=ALAMF*TORAD
      ALAMR=ALAMR*TORAD
      ROTANG=ROTANG*TORAD
      DO 253 J=1,JMAX
      THETA(J)=THETA(J)*TORAD
      THETOB(J) = THETOB(J)*TORAD
      OBANGS(J)=OBANGS(J)*TORAD
  253 CONTINUE
C
C        IF SHARP LEADING-EDGE IS USED, SET JLE1A AND JLE2A
      IF(NLETYP.EQ.1)  GO TO 260
      WRITE(6,300)
      GO TO (266,267,268),NTETYP
  266 NSKP=4
      JST=3
      GO TO 269
  267 NSKP=3
      JST=2
      GO TO 269
  268 NSKP=3
      JST=4
C
  269 DO 262 L=1,2
      JTETM=JTETOP-NSKP
C
      DO 261 J=JST,JTETM,NSKP
      J2=J+NSKP
      IF(AIRFY(J2).GT.0.)  GO TO 263
  261 CONTINUE
      WRITE(6,381)
      FATAL=.TRUE.
      GO TO 260
C
  263 JLE1A(L)=J
      JLE2A(L)=J2
      IF(NLETYP.EQ.2)  WRITE(6,383)  L,J
      IF(NLETYP.EQ.3)  WRITE(6,382)  L,J,L,J2
      NSKP=1
      JST=1
      IF(NTETYP.EQ.3)  JST=2
  262 CONTINUE
  260 CONTINUE
C
C        IF ANY FATAL ERRORS HAVE BEEN FOUND IN THE INPUT, STOP.
      IF(FATAL)  WRITE(6,375)
      IF(FATAL)  STOP
      RETURN
C
  100 FORMAT(51H0FATAL ERROR.  VALUE FOR  JAIRF IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  101 FORMAT(51H0FATAL ERROR.  VALUE FOR  JMAX  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  102 FORMAT(51H0FATAL ERROR.  VALUE FOR  KMAX  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  103 FORMAT(51H0FATAL ERROR.  VALUE FOR NTETYP IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  104 FORMAT(51H0FATAL ERROR.  VALUE FOR NAIRF  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  105 FORMAT(51H0FATAL ERROR.  VALUE FOR NIBDST IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  106 FORMAT(51H0FATAL ERROR.  VALUE FOR   TR   IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  108 FORMAT(51H0FATAL ERROR.  BAD INPUT VALUES FOR XLE AND/OR XTE.,
     1  28H  XLE MUST BE LESS THAN XTE.)
  109 FORMAT(51H0FATAL ERROR.  VALUE FOR NOBSHP IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  110 FORMAT(51H0FATAL ERROR.  VALUE FOR NOBDST IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  111 FORMAT(47H0FATAL ERROR.  INCONSISTANCY IN THE VALUES FOR ,
     1  28HRADOB, XOBCNT, XLE, AND XTE./29H THE AIRFOIL MUST FIT INSIDE
     *,
     2  22HOF THE OUTER BOUNDARY.)
  112 FORMAT(53H0FATAL ERROR.  BAD INPUT VALUE FOR XLEFT.  XLEFT MUST,
     1  18H BE LESS THAN XLE.)
  113 FORMAT(54H0FATAL ERROR. BAD INPUT VALUE FOR XRIGHT.  XRIGHT MUST,
     1  21H BE GREATER THAN XTE.)
  114 FORMAT(51H0FATAL ERROR.  VALUE FOR YBOTOM IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  115 FORMAT(51H0FATAL ERROR.  VALUE FOR  YTOP  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  116 FORMAT(40H0FATAL ERROR. BAD INPUT VALUE FOR RCORN./
     1  54H RCORN MUST BE GREATER THAN ZERO, AND MUST MAKE SENSE ,
     2  41HWITH RESPECT TO THE SIZE OF THE RECTANGLE/12H DEFINED BY ,
     3  32HXLEFT, XRIGHT, YBOTOM, AND YTOP.)
  117 FORMAT(51H0FATAL ERROR.  VALUE FOR ALAMF  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  118 FORMAT(51H0FATAL ERROR.  VALUE FOR ALAMR  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  119 FORMAT(53H0FATAL ERROR.  VALUE FOR NORDA(1) IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  120 FORMAT(53H0FATAL ERROR.  VALUE FOR NORDA(2) IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  121 FORMAT(54H0FATAL ERROR.  VALUE FOR MAXITA(1) IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  122 FORMAT(54H0FATAL ERROR.  VALUE FOR MAXITA(1) IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  124 FORMAT(51H0FATAL ERROR.  VALUE FOR  NPLT  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  125 FORMAT(41H0FATAL ERROR.  BAD INPUT VALUE FOR JAIRF./10H IF NAIRF ,
     1  51HEQUALS 5, JAIRF MUST BE GREATER THAN OR EQUAL TO 4.)
  126 FORMAT(41H0FATAL ERROR.  BAD INPUT VALUE FOR NAIRF./6H WITH ,
     1  43HNIBDST EQUAL TO 5, NAIRF MUST ALSO EQUAL 5.)
  127 FORMAT(53H0SINCE NTETYP INDICATES THAT THIS IS AN O-TYPE GRID, ,
     1  34HJTEBOT AND JTETOP HAVE BEEN SET TO/22H THEIR DEFAULT VALUES ,
     2  27HOF 1 AND JMAX, RESPECTIVLY.)
  128 FORMAT(50H0FATAL ERROR.  BAD INPUT VALUES FOR JTEBOT AND/OR ,
     1  7HJTETOP.)
  129 FORMAT(53H0JMAX IS NOT A MULTIPLE OF 4.  THUS, FOR NTETYP = 1, ,
     1  26HSEQUENCING CANNOT BE USED.)
  130 FORMAT(53H0JMAX IS NOT A MULTIPLE OF 3.  THUS, FOR NTETYP = 2, ,
     1  26HSEQUENCING CANNOT BE USED.)
  131 FORMAT(55H0JMAX IS NOT A MULTIPLE OF 3 PLUS 1.  THUS, FOR NTETYP ,
     1  31H= 3, SEQUENCING CANNOT BE USED.)
  132 FORMAT(46H0GRID SEQUENCING WILL BE USED.  THERE WILL BE ,
     1  46HA COARSE SOLUTION FOLLOWED BY A FINE SOLUTION.)
  133 FORMAT(52H0FATAL ERROR.  VALUE FOR   JD    IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  137 FORMAT(49H0FATAL ERROR.  FOR CASCADE GRIDS (NOBSHP=3 OR 5) ,
     1  40HALAMF AND ALAMR MAY NOT BOTH EQUAL ZERO.)
  139 FORMAT(50H0FATAL ERROR.  VALUE FOR WAKEP IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  140 FORMAT(50H0FATAL ERROR.  BAD INPUT VALUES FOR NTETYP AND/OR ,
     1  7HNOBSHP./49H IF NTETYP INDICATES AN O-TYPE GRID, NOBSHP MUST ,
     2  17HEQUAL 1, 2, OR 3./36H IF NTETYP INDICATES A C-TYPE GRID, ,
     3  25HNOBSHP MUST EQUAL 4 OR 5.)
  141 FORMAT(51H0FATAL ERROR.  VALUE FOR NOBCAS IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  142 FORMAT(50H0FATAL ERROR.  BAD INPUT VALUES FOR NTETYP AND/OR ,
     1  7HNOBDST./49H FOR CASCADE OUTER BOUNDARIES, NOBDST MUST EQUAL ,
     2  2H1.)
  143 FORMAT(41H0FATAL ERROR.  BAD INPUT VALUE FOR ALAMR./5H FOR ,
     1  57HNOBSHP INDICATING THAT THIS IS NOT A CASCADE GRID, ALAMR ,
     2  16HMUST EQUAL ZERO.)
  144 FORMAT(52H0SINCE NIBDST=5, XLE AND XTE HAVE BEEN REDEFINED TO ,
     1  46HCORRESPOND TO GIVEN AIRFOIL COORDINATES.  XLE=,F13.6,5X,
     2  4HXTE=,F13.6)
  146 FORMAT(51H0FATAL ERROR.  VALUE FOR NLETYP IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  152 FORMAT(51H0FATAL ERROR.  VALUE FOR  BINN  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  153 FORMAT(51H0FATAL ERROR.  VALUE FOR JCAMBR IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  154 FORMAT(51H0FATAL ERROR.  VALUE FOR XTFRAC IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  155 FORMAT(51H0FATAL ERROR.  VALUE FOR ROTANG IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  158 FORMAT(51H0FATAL ERROR.  VALUE FOR OMEGA  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  159 FORMAT(51H0FATAL ERROR.  VALUE FOR OMEGP  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  160 FORMAT(51H0FATAL ERROR.  VALUE FOR OMEGQ  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  161 FORMAT(51H0FATAL ERROR.  VALUE FOR OMEGR  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  162 FORMAT(51H0FATAL ERROR.  VALUE FOR OMEGS  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  163 FORMAT(51H0FATAL ERROR.  VALUE FOR  PLIM  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  164 FORMAT(51H0FATAL ERROR.  VALUE FOR  QLIM  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  165 FORMAT(51H0FATAL ERROR.  VALUE FOR  RLIM  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  166 FORMAT(51H0FATAL ERROR.  VALUE FOR  SLIM  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  167 FORMAT(51H0FATAL ERROR.  VALUE FOR   NDS  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  168 FORMAT(51H0FATAL ERROR.  VALUE FOR   DSI  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  169 FORMAT(51H0FATAL ERROR.  VALUE FOR DSOBI  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  170 FORMAT(51H0FATAL ERROR.  VALUE FOR THETAI IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  171 FORMAT(51H0FATAL ERROR.  VALUE FOR THOBI  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  172 FORMAT(51H0FATAL ERROR.  VALUE FOR  AAAI  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  173 FORMAT(51H0FATAL ERROR.  VALUE FOR  BBBI  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  174 FORMAT(51H0FATAL ERROR.  VALUE FOR  CCCI  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  175 FORMAT(51H0FATAL ERROR.  VALUE FOR  DDDI  IS OUTSIDE OF ALLOW,
     1  11HABLE RANGE.)
  176 FORMAT(55H FATAL ERROR.  VALUE FOR TEOPEN IS OUTSIDE OF ALLOWABLE,
     1  7H RANGE.)
  193 FORMAT(53H0FATAL ERROR.  VALUE FOR NOUT IS OUTSIDE OF ALLOWABLE,
     1  7H RANGE.)
  194 FORMAT(49H0SEQUENCING WILL NOT BE USED BECAUSE MAXITA(1)=0.)
  195 FORMAT(51H0SEQUENCING WILL NOT BE USED BECAUSE KMAX IS NOT A ,
     1  21HMULTIPLE OF 3 PLUS 1.)
  300 FORMAT(1H )
  301 FORMAT(49H WARNING.  INPUT ARRAY   DS   MAY NOT BE SMOOTH. ,
     1  14H CHECK     DS(,I3,2H)=,F13.6)
  302 FORMAT(49H WARNING.  INPUT ARRAY  DSOB  MAY NOT BE SMOOTH. ,
     1  14H CHECK   DSOB(,I3,2H)=,F13.6)
  303 FORMAT(49H WARNING.  INPUT ARRAY THETA  MAY NOT BE SMOOTH. ,
     1  14H CHECK  THETA(,I3,2H)=,F13.6)
  304 FORMAT(49H WARNING.  INPUT ARRAY THETOB MAY NOT BE SMOOTH. ,
     1  14H CHECK THETOB(,I3,2H)=,F13.6)
  305 FORMAT(49H WARNING.  INPUT ARRAY  DIST  MAY NOT BE SMOOTH. ,
     1  14H CHECK   DIST(,I3,2H)=,F13.6)
  306 FORMAT(49H WARNING.  INPUT ARRAY AIRFX  MAY NOT BE SMOOTH. ,
     1  14H CHECK  AIRFX(,I3,2H)=,F13.6)
  307 FORMAT(49H WARNING.  INPUT ARRAY AIRFY  MAY NOT BE SMOOTH. ,
     1  14H CHECK  AIRFY(,I3,2H)=,F13.6)
  308 FORMAT(49H WARNING.  INPUT ARRAY CAMBRX MAY NOT BE SMOOTH. ,
     1  14H CHECK CAMBRX(,I3,2H)=,F13.6)
  309 FORMAT(49H WARNING.  INPUT ARRAY CAMBRY MAY NOT BE SMOOTH. ,
     1  14H CHECK CAMBRY(,I3,2H)=,F13.6)
  310 FORMAT(49H WARNING.  INPUT ARRAY OBANGS MAY NOT BE SMOOTH. ,
     1  14H CHECK OBANGS(,I3,2H)=,F13.6)
  311 FORMAT(49H WARNING.  INPUT ARRAY  AAA   MAY NOT BE SMOOTH. ,
     1  14H CHECK    AAA(,I3,2H)=,F13.6)
  312 FORMAT(49H WARNING.  INPUT ARRAY  BBB   MAY NOT BE SMOOTH. ,
     1  14H CHECK    BBB(,I3,2H)=,F13.6)
  313 FORMAT(49H WARNING.  INPUT ARRAY  CCC   MAY NOT BE SMOOTH. ,
     1  14H CHECK    CCC(,I3,2H)=,F13.6)
  314 FORMAT(49H WARNING.  INPUT ARRAY  DDD   MAY NOT BE SMOOTH. ,
     1  14H CHECK    DDD(,I3,2H)=,F13.6)
  316 FORMAT(41H WARNING.  CHECK XMIN, XMAX, YMIN, YMAX =,4F13.6,
     1  8H  FOR N=,I3)
  321 FORMAT(32H FATAL ERROR.  VALUE FOR     DS(,I3,2H)=,F13.6,5X,
     1  26HIS OUT OF ALLOWABLE RANGE.)
  322 FORMAT(32H FATAL ERROR.  VALUE FOR   DSOB(,I3,2H)=,F13.6,5X,
     1  26HIS OUT OF ALLOWABLE RANGE.)
  323 FORMAT(32H FATAL ERROR.  VALUE FOR  THETA(,I3,2H)=,F13.6,5X,
     1  26HIS OUT OF ALLOWABLE RANGE.)
  324 FORMAT(32H FATAL ERROR.  VALUE FOR THETOB(,I3,2H)=,F13.6,5X,
     1  26HIS OUT OF ALLOWABLE RANGE.)
  325 FORMAT(32H FATAL ERROR.  VALUE FOR   DIST(,I3,2H)=,F13.6,5X,
     1  26HIS OUT OF ALLOWABLE RANGE.)
  328 FORMAT(32H FATAL ERROR.  VALUE FOR CAMBRX(,I3,2H)=,F13.6,5X,
     1  26HIS OUT OF ALLOWABLE RANGE.)
  330 FORMAT(32H FATAL ERROR.  VALUE FOR OBANGS(,I3,2H)=,F13.6,5X,
     1  26HIS OUT OF ALLOWABLE RANGE.)
  331 FORMAT(32H FATAL ERROR.  VALUE FOR    AAA(,I3,2H)=,F13.6,5X,
     1  26HIS OUT OF ALLOWABLE RANGE.)
  332 FORMAT(32H FATAL ERROR.  VALUE FOR    BBB(,I3,2H)=,F13.6,5X,
     1  26HIS OUT OF ALLOWABLE RANGE.)
  333 FORMAT(32H FATAL ERROR.  VALUE FOR    CCC(,I3,2H)=,F13.6,5X,
     1  26HIS OUT OF ALLOWABLE RANGE.)
  334 FORMAT(32H FATAL ERROR.  VALUE FOR    DDD(,I3,2H)=,F13.6,5X,
     1  26HIS OUT OF ALLOWABLE RANGE.)
  340 FORMAT(55H FATAL ERROR.  ARRAY   DIST IS NOT MONOTONIC.  ELEMENT ,
     1  I3,34H FAILS TO BE GREATER THAN ELEMENT ,I3,1H.)
  341 FORMAT(55H FATAL ERROR.  ARRAY CAMBRX IS NOT MONOTONIC.  ELEMENT ,
     1  I3,34H FAILS TO BE GREATER THAN ELEMENT ,I3,1H.)
  342 FORMAT(55H FATAL ERROR.  ARRAY OBANGS IS NOT MONOTONIC.  ELEMENT ,
     1  I3,34H FAILS TO BE GREATER THAN ELEMENT ,I3,1H.)
  343 FORMAT(53H FOR THIS PARTICULAR CHOICE OF NOBSHP, XRIGHT, YTOP, ,
     1  11HAND YBOTOM,/45H THE ELEMENTS OF ARRAY OBANGS MUST BE BETWEEN,
     2  F9.2,4H AND,F9.2,9H DEGREES.)
  344 FORMAT(55H0FATAL ERROR.  THE ELEMENTS OF ARRAY AIRFX ARE ALL THE ,
     1  5HSAME.)
  353 FORMAT(50H0FATAL ERROR. FOR SHARP LEADING-EDGE, NIBDST MUST ,
     1  8HEQUAL 5.)
  373 FORMAT(47H WARNING.  INPUT ARRAY XOB MAY NOT BE SMOOTH.  ,
     1  11H CHECK XOB(,I3,2H)=,F13.6)
  374 FORMAT(47H WARNING.  INPUT ARRAY YOB MAY NOT BE SMOOTH.  ,
     1  11H CHECK YOB(,I3,2H)=,F13.6)
  375 FORMAT(48H0FATAL ERROR(S) FOUND IN INPUT VARIABLES.  EXIT.)
  376 FORMAT(1H1,/////////  20X,101(1HX)/ 20X,1HX,99X,1HX/
     1  20X,1HX,99X,1HX/  20X,1HX,99X,1HX/
     2  20X,1HX,T49,41HAPPLIED COMPUTATIONAL AERODYNAMICS BRANCH,
     3  T121,1HX/ 20X,1HX,99X,1HX/  20X,1HX,99X,1HX/
     4  20X,1HX,T55,27HNASA - AMES RESEARCH CENTER,T121,1HX/
     5  20X,1HX,99X,1HX/  20X,1HX,99X,1HX/
     6  20X,1HX,T57,23HMOFFETT FIELD, CA 94035,T121,1HX/
     7  20X,1HX,99X,1HX/  20X,1HX,99X,1HX/  20X,1HX,T47,
     8  47H GGGGG    RRRRRR       A      PPPPPP    EEEEEEE,
     9  T121,1HX/ 20X,1HX,T47,
     A  47HGG        RR   RR     A A     PP   PP   EE     ,
     B  T121,1HX/ 20X,1HX,T47,
     C  47HGG GGGG   RRRRRR     A   A    PPPPPP    EEEEE  ,
     D  T121,1HX/ 20X,1HX,T47,
     E  47HGG G GG   RR RR     AAAAAAA   PP        EE     ,
     F  T121,1HX/ 20X,1HX,T47,
     G  47HGG   GG   RR  RR    AA   AA   PP        EE     ,
     H  T121,1HX)
  377 FORMAT(20X,1HX,T47,
     I  47H GGGG G   RR   RR   AA   AA   PP        EEEEEEE,
     J  T121,1HX/
     Y  20X,1HX,99X,1HX/  20X,1HX,T47,28H(GRIDS ABOUT AIRFOILS USING ,
     Z  17HPOISSON EQUATION),T121,1HX/
     K  20X,1HX,99X,1HX/  20X,1HX,99X,1HX/  20X,1HX,99X,1HX/
     L  20X,1HX,T31,39HTHIS PROGRAM COMPUTES TWO-DIMENSIONAL, ,
     M  36HCURVALINEAR, FINITE DIFFERENCE GRIDS,T121,1HX/
     O  20X,1HX,T31,37HABOUT AIRFOILS BY SOLVING THE POISSON,
     P  42H DIFFERENTIAL EQUATION TO ACHIEVE CONTROL ,T121,1HX/
     R 20X,1HX,T50,41HOF CLUSTERING AND ANGLES NEAR BOUNDARIES.,T121,1HX
     */
     S  20X,1HX,99X,1HX/  20X,1HX,99X,1HX/  20X,1HX,99X,1HX/
     T  20X,1HX,T65,10HWRITTEN BY,T121,1HX/
     U  20X,1HX,99X,1HX/  20X,1HX,99X,1HX/
     V  20X,1HX,T62,17HREESE L. SORENSON,T121,1HX/
     W  20X,1HX,99X,1HX/  20X,1HX,99X,1HX/  20X,1HX,99X,1HX/
     X  20X,101(1HX)/////)
  378 FORMAT(43H MAXITA(1)  ADDED TO MAXITA(2).  MAXITA(2)=,I5/
     1  40H NORDA(1)  ADDED TO NORDA(2).  NORDA(2)=,I2)
  379 FORMAT(54H0FATAL ERROR.  FOR NIBDST LESS THAN 5, THUS INDICATING,
     1  51H THAT INTERPOLATION IS TO BE DONE, THE AIRFOIL MUST/
     2  52H GO ALL THE WAY AROUND FROM THE TRAILING EDGE TO THE,
     3  49H TRAILING EDGE AGAIN.  THIS CONDITION IS NOT MET.)
  380 FORMAT(51H0FATAL ERROR.  FOR NLETYP GREATER THAN 1, JMAX MUST,
     1  13H EQUAL JAIRF.)
  381 FORMAT(54H FATAL ERROR.  FAILURE TO FIND LEADING-EDGE OF AIRFOIL,
     1  19H FOR NLETYP=3 OR 3.)
  382 FORMAT(7H JLE1A(,I1,2H)=,I3,10X,6HJLE2A(,I1,2H)=,I3)
  383 FORMAT(7H JLE1A(,I1,2H)=,I3)
      END
*DECK INNER
      SUBROUTINE INNER
C
C        THIS SUBROUTINE DISTRIBUTES POINTS ON THE INNER (AIRFOIL)
C        BOUNDARY.
C
*CALL GRIDA
C
*CALL GRIDB
C
*CALL GRIDC
C
      DIMENSION R(241),ADIST(170),S(170),BDIST(241),GKX(161),GKY(161),
     1  D64X(51),D64Y(51),DCPM(241),YC(170),YCP(170)
      DIMENSION GKY1(72),GKY2(89),DCPM1(60),DCPM2(60),DCPM3(60),
     1  DCPM4(61)
C
      EQUIVALENCE (GKY(1),GKY1(1)),(GKY(73),GKY2(1)),(DCPM(1),DCPM1(1)),
     1  (DCPM(61),DCPM2(1)),(DCPM(121),DCPM3(1)),(DCPM(181),DCPM4(1))
C
      DATA GKX/
     1 1.     ,.99953,.99810,.99572,.99237,.98806,.98278,.97654,.96934,
     2  .96120,.95212,.94212,.93122,.91943,.90679,.89331,.87902,.86398,
     3  .84821,.83177,.81471,.79713,.77909,.76068,.74198,.72309,.70405,
     4  .68491,.66571,.64646,.62720,.60792,.58864,.56937,.55012,.53091,
     5  .51174,.49263,.47360,.45467,.43584,.41714,.39859,.38020,.36200,
     6  .34401,.32624,.30873,.29149,.27455,.25793,.24165,.22574,.21022,
     7  .19512,.18044,.16622,.15247,.13920,.12644,.11420,.10250,.09134,
     8  .08075,.07074,.06131,.05250,.04431,.03676,.02987,.02366,.01814,
     9  .01336,.00932,.00605,.00351,.00165,.00049,.00000,.00003,.00056,
     A  .00171,.00365,.00640,.00993,.01415,.01906,.02461,.03080,.03761,
     B  .04503,.05304,.06163,.07080,.08057,.09101,.10215,.11398,.12642,
     C  .13946,.15306,.16720,.18183,.19694,.21251,.22850,.24488,.26164,
     D  .27873,.29615,.31385,.33181,.35000,.36839,.38695,.40565,.42446,
     E  .44335,.46229,.48125,.50019,.51910,.53793,.55666,.57526,.59371,
     F  .61198,.63006,.64794,.66559,.68302,.70019,.71710,.73373,.75007,
     G  .76609,.78179,.79714,.81214,.82675,.84096,.85475,.86810,.88099,
     H  .89338,.90526,.91660,.92735,.93749,.94699,.95580,.96390,.97125,
     I  .97783,.98360,.98854,.99262,.99582,.99813,.99953,1./
C
      DATA GKY1/
     1   0.    , .00006, .00022, .00047, .00080, .00119, .00163, .00208,
     2   .00255, .00300, .00342, .00378, .00407, .00427, .00435, .00428,
     3   .00405, .00362, .00296, .00206, .00088,-.00059,-.00236,-.00442,
     4  -.00675,-.00930,-.01202,-.01485,-.01774,-.02065,-.02352,-.02634,
     5  -.02907,-.03168,-.03418,-.03653,-.03873,-.04077,-.04265,-.04436,
     6  -.04589,-.04726,-.04844,-.04946,-.05053,-.05097,-.05147,-.05181,
     7  -.05198,-.05200,-.05186,-.05156,-.05112,-.05052,-.04978,-.04889,
     8  -.04786,-.04669,-.04538,-.04394,-.04237,-.04067,-.03886,-.03693,
     9  -.03491,-.03280,-.03061,-.02834,-.02602,-.02364,-.02122,-.01877/
      DATA GKY2/
     A  -.01626,-.01370,-.01105,-.00827,-.00538,-.00240, .00074, .00399,
     B   .00717, .01015, .01296, .01573, .01852, .02133, .02415, .02696,
     C   .02973, .03245, .03510, .03766, .04009, .04237, .04444, .04629,
     D   .04800, .04963, .05118, .05266, .05408, .05543, .05670, .05791,
     F   .05903, .06007, .06104, .06192, .06271, .06341, .06403, .06455,
     G   .06499, .06532, .06556, .06571, .06575, .06569, .06553, .06526,
     H   .06489, .06439, .06379, .06305, .06219, .06118, .06003, .05873,
     I   .05727, .05566, .05390, .05200, .04996, .04780, .04552, .04315,
     J   .04069, .03816, .03558, .03296, .03032, .02769, .02508, .02252,
     K   .02003, .01762, .01533, .01316, .01114, .00928, .00759, .00607,
     L   .00472, .00355, .00256, .00175, .00109, .00060, .00026, .00006,
     M   0./
C
      DATA D64X/
     1  1.    ,.94947,.89896,.84852,.79849,.74874,.69892,.64915,.59943,
     2  .54975,.50011,.45050,.40090,.35129,.30166,.25200,.20230,.15252,
     3  .10263,.07770,.05251,.02724,.01441,.00918,.00650,.0    ,.00350,
     4  .00582,.01059,.02276,.04749,.07230,.09737,.14748,.19770,.24800,
     5  .29834,.34871,.39910,.44950,.49989,.55025,.60057,.65085,.70108,
     6  .75126,.80151,.85148,.90104,.95053,1./
C
      DATA D64Y/
     1  -.00021,-.00048,-.00076,-.00132,-.00229,-.00460,-.00760,-.01086,
     2  -.01418,-.01736,-.02024,-.02266,-.02436,-.02518,-.02537,-.02499,
     3  -.02406,-.02244,-.01996,-.01919,-.01592,-.01251,-.00969,-.00796,
     4  -.00678, .00000, .00902, .01112, .01451, .02095, .03034, .03865,
     5   .04380, .05366, .06126, .06705, .07131, .07414, .07552, .07522,
     6   .07344, .07040, .06624, .06106, .05490, .04780, .03967, .03018,
     7   .02038, .01028, .00021/
C
      DATA DCPM1/
     1.000000000,.000293863,.000692701,.001244258,.001966212,.002865780,
     2.003946105,.005208424,.006652930,.008279162,.010086195,.012072753,
     3.014237271,.016577934,.019092707,.021779350,.024635440,.027658375,
     4.030845390,.034193561,.037699818,.041360949,.045173606,.049134314,
     5.053239480,.057485392,.061868234,.066384086,.071028934,.075798675,
     6.080689121,.085696010,.090815009,.096041719,.101371685,.106800398,
     7.112323303,.117935803,.123633267,.129411033,.135264418,.141188715,
     8.147179209,.153231172,.159339874,.165500585,.171708584,.177959155,
     9.184247601,.190569242,.196919419,.203293502,.209686891,.216095018,
     A.222513353,.228937407,.235362733,.241784933,.248199654,.254602598/
      DATA DCPM2/
     B.260989519,.267356230,.273698599,.280012558,.286294099,.292539281,
     C.298744226,.304905126,.311018242,.317079904,.323086515,.329034550,
     D.334920558,.340741164,.346493070,.352173051,.357777963,.363304739,
     E.368750394,.374112020,.379386791,.384571964,.389664875,.394662949,
     F.399563689,.404364687,.409063621,.413658255,.418146442,.422526127,
     G.426795345,.430952227,.434994998,.438921985,.442731616,.446422428,
     H.449993066,.453442297,.456769009,.459972223,.463051105,.466004974,
     I.468833322,.471535831,.474112397,.476563158,.478888537,.481089285,
     J.483166543,.485121923,.486957607,.488676479,.490282289,.491779869,
     K.493175398,.494476716,.495693675,.496838478,.497925903,.498973269/
      DATA DCPM3/
     L.500000000,.501026731,.502074097,.503161522,.504306325,.505523284,
     M.506824602,.508220131,.509717711,.511323521,.513042393,.514878077,
     N.516833457,.518910715,.521111463,.523436842,.525887603,.528464169,
     O.531166678,.533995026,.536948895,.540027777,.543230991,.546557703,
     P.550006934,.553577572,.557268384,.561078015,.565005002,.569047773,
     Q.573204655,.577473873,.581853558,.586341745,.590936379,.595635313,
     R.600436311,.605337051,.610335125,.615428036,.620613209,.625887980,
     S.631249606,.636695261,.642222037,.647826949,.653506930,.659258836,
     T.665079442,.670965450,.676913485,.682920096,.688981758,.695094874,
     U.701255774,.707460719,.713705901,.719987442,.726301401,.732643770/
      DATA DCPM4/
     V.739010481,.745397402,.751800346,.758215067,.764637267,.771062593,
     W.777486647,.783904982,.790313109,.796706498,.803080581,.809430758,
     X.815752399,.822040845,.828291416,.834499415,.840660126,.846768828,
     Y.852820791,.858811285,.864735582,.870588967,.876366733,.882064197,
     Z.887676697,.893199602,.898628315,.903958281,.909184991,.914303990,
     1.919310879,.924201325,.928971066,.933615914,.938131766,.942514608,
     2.946760520,.950865686,.954826394,.958639051,.962300182,.965806439,
     3.969154610,.972341625,.975364560,.978220650,.980907293,.983422066,
     4.985762729,.987927247,.989913805,.991720838,.993347070,.994791576,
     5.996053895,.997134220,.998033788,.998755742,.999307299,.999706137,
     6  1./
C
C        *************************************************************
C        *                                                           *
C        *     PART 1:  SPECIAL TREATMENT FOR THE CASE (NIBDST=5)    *
C        *     WHERIN AIRFOIL CARDS HAVE BEEN READ AND THEY ARE TO   *
C        *     BE USED EXACTLY AS READ.                              *
C        *                                                           *
C        *************************************************************
      IF(NIBDST.NE.5)  GO TO 3
      L=0
      DO 1 J=JTEBOT,JTETOP
      L=L+1
      X(J,1)=AIRFX(L)
    1 Y(J,1)=AIRFY(L)
      GO TO 37
C
C        *************************************************************
C        *                                                           *
C        *     PART 2:  FIND ADIST, THE DISTRIBUTION FUNCTION WE ARE *
C        *     GOING TO USE.                                         *
C        *                                                           *
C        *************************************************************
    3 GO TO (4,6,7,5),NIBDST
C
C        USE PRE-STORED CIRCLE-PLANE MAPPING DISTRIBUTION.
    4 JDIST=241
      DO 8 J=1,JDIST
    8 DIST(J)=DCPM(J)
      GO TO 5
C
C        CALCULATE DISTRIBUTION FROM FORMULA.
    6 JDIST=241
      BIP=BINN+1.
      BIM=BINN-1.
      C1=ALOG(BIP/BIM)
      BIP3=BINN+3.
      BIM3=BINN-3.
      DT1=1./FLOAT(JDIST-1)
      T1=-DT1
C
      DO 15 J=1,JDIST
      T1=T1+DT1
      IF(T1.GT..5)  GO TO 16
      EXPA=EXP(C1*(4.*T1-1.))
      DIST(J)=(BIM-BIP*EXPA)/(-4.*(EXPA+1.))
      GO TO 15
   16 EXPA=EXP(C1*(4.*T1-3.))
      DIST(J)=(BIM3-BIP3*EXPA)/(-4.*(EXPA+1.))
   15 CONTINUE
      GO TO 5
C
C        TAKE DISTRIBUTOIN FROM AIRFOIL DATA CARDS.
    7 JDIST=JAIRF
      DIST(1)=0.
      DO 9 J=2,JDIST
    9 DIST(J)=DIST(J-1)+SQRT((AIRFX(J)-AIRFX(J-1))**2+(AIRFY(J)-
     1  AIRFY(J-1))**2)
      RLAST=1./DIST(JDIST)
      DO 10 J=2,JDIST
   10 DIST(J)=DIST(J)*RLAST
C
C        IF NIBDST=4 USE USE DISTRIBUTION FUNCTION AS SUPPLIED BY
C        USER ON CARDS.
    5 DR=1./FLOAT(JDIST-1)
      R(1)=0.
      DO 11 J=2,JDIST
   11 R(J)=R(J-1)+DR
C
C        DICIDE WHERE THE DISTRIBUTION SHOULD START AND END, DEPENDING
C        ON WHETHER OR NOT THERE ARE POINTS AT THE TRAILING EDGE.
      GO TO (12,13,14),NTETYP
   12 JBOD=JMAX
      SST=0.
      SND=FLOAT(JBOD-1)
      GO TO 17
   13 JBOD=JMAX
      SST=.5
      SND=FLOAT(JBOD)-.5
      GO TO 17
   14 JBOD=JTETOP-JTEBOT+1
      SST=0.
      SND=FLOAT(JBOD)
   17 DELS=(SND-SST)/FLOAT((JBOD-1)*JBOD)
      S(1)=SST*DELS
      DO 18 J=2,JBOD
   18 S(J)=S(J-1)+DELS
C
      CALL CSPLIN(S,ADIST,DUM,R,DIST,1,JBOD,1,JDIST,0)
C
C        *************************************************************
C        *                                                           *
C        *     PART3:  FIND AIRFX AND AIRFY, THE AIRFOIL SHAPE WE    *
C        *     ARE GOING TO USE.                                     *
C        *                                                           *
C        *************************************************************
      GO TO (24,22,19,26,28),NAIRF
C
C        USE TRUE NACA 00XX AIRFOIL, WHICH HAS OPEN TRAILING-EDGE.
   24 XT=1.
      GO TO 23
C
C        USE NACA 00XX AIRFOIL, MODIFIED BY EXTRAPOLATING THE ANALYTIC
C        DEFINING FUNCTION TO PROVIDE A CLOSED TRAILING-EDGE.
   22 XT=1.008930411365
C
   23 RXT=1./XT
      YSIGN=-1.
      JAIRF=241
      DELZ=1./FLOAT(JAIRF-1)
      Z=-DELZ
      XX=1.5
C
      DO 25 J=1,JAIRF
      Z=Z+DELZ
      XL=XX
      XX=XT*((.5*(1.+COS(TWOPI*Z)))**1.5)
      IF(XX.GE.XL)  YSIGN=1.
      AIRFX(J)=XX*RXT
   25 AIRFY(J)=SIGN(5.*TR*(.2969*SQRT(XX)-.126*XX-.3516*XX**2+.2843*XX
     1  **3-.1015*XX**4),YSIGN)
      GO TO 20
C
C        USE GARABEDIAN-KORN AIRFOIL.
   19 JAIRF=161
      DO 21 J=1,JAIRF
      AIRFX(J)=GKX(J)
   21 AIRFY(J)=GKY(J)
      GO TO 20
C
C        USE NACA 64A410 AIRFOIL.
   26 JAIRF=51
      DO 27 J=1,JAIRF
      AIRFX(J)=D64X(J)
   27 AIRFY(J)=D64Y(J)
      GO TO 20
C
C        USE AIRFOIL SHAPE AS READ FROM DATA CARDS.  NORMALIZE IT TO GO
C        FROM 0 TO 1 IF IT DOES NOT ALREADY.
   28 XMX=AIRFX(1)
      XMN=XMX
      DO 29 J=2,JAIRF
      IF(AIRFX(J).GT.XMX)  XMX=AIRFX(J)
      IF(AIRFX(J).LT.XMN)  XMN=AIRFX(J)
   29 CONTINUE
C
      IF(XMN.EQ.0..AND.XMX.EQ.1.)  GO TO 20
      XRAT=1./(XMX-XMN)
      DO 30 J=1,JAIRF
      AIRFX(J)=(AIRFX(J)-XMN)*XRAT
   30 AIRFY(J)=AIRFY(J)*XRAT
C
C        FIND BDIST, THE ARCLENGTH DISTRIBUTION OF THE GIVEN AIRFX AND
C        AIRFY.
   20 BDIST(1)=0.
      DO 31 J=2,JAIRF
   31 BDIST(J)=BDIST(J-1)+SQRT((AIRFX(J)-AIRFX(J-1))**2+(AIRFY(J)
     1  -AIRFY(J-1))**2)
      RLAST=1./BDIST(JAIRF)
      DO 43 J=2,JAIRF
   43 BDIST(J)=BDIST(J)*RLAST
      TEOPEN=(AIRFY(JAIRF)-AIRFY(1))*(XTE-XLE)
C
C        *************************************************************
C        *                                                           *
C        *     PART4:  GIVEN AIRFX AND AIRFY AS FUNCTIONS OF BDIST,  *
C        *     FIND THE ACTUAL FOUNDARY POINTS AS FUNCTINNS OF ADIST.*
C        *                                                           *
C        *************************************************************
      CALL CSPLIN(ADIST,X(JTEBOT,1),DUM,BDIST,AIRFX,1,JBOD,1,JAIRF,0)
C
      CALL CSPLIN(ADIST,Y(JTEBOT,1),DUM,BDIST,AIRFY,1,JBOD,1,JAIRF,0)
C
C        *************************************************************
C        *                                                           *
C        *     PART 5:  IF CALLED-FOR, CAMBER THE AIRFOIL USING A    *
C        *     USER-SUPPLIED CAMBER-CURVE.                           *
C        *                                                           *
C        *************************************************************
      IF(JCAMBR.EQ.0)  GO TO 32
      CALL CSPLIN(X(JTEBOT,1),YC,YCP,CAMBRX,CAMBRY,1,JBOD,1,JCAMBR,1)
      L=0
      XT=1.5
      DO 33 J=JTEBOT,JTETOP
      L=L+1
      TH=ATAN(YCP(L))
      YT=ABS(Y(J,1))
      XTL=XT
      XT=X(J,1)
      IF(Y(J,1).GT.0)  GO TO 34
      X(J,1)=XT+YT*SIN(TH)
      Y(J,1)=YC(L)-YT*COS(TH)
      GO TO 33
   34 X(J,1)=XT-YT*SIN(TH)
      Y(J,1)=YC(L)+YT*COS(TH)
   33 CONTINUE
C
C        *************************************************************
C        *                                                           *
C        *     PART 6:  DO A LINEAR TRANSFORMATION ON THE AIRFOIL    *
C        *     TO MAKE IT TO GO FROM XLE TO XTE.                     *
C        *                                                           *
C        *************************************************************
   32 IF(XLE.EQ.0..AND.XTE.EQ.1.)  GO TO 35
      XRATIO=XTE-XLE
      DO 36 J=JTEBOT,JTETOP
      X(J,1)=XLE+X(J,1)*XRATIO
   36 Y(J,1)=Y(J,1)*XRATIO
C
C        *************************************************************
C        *                                                           *
C        *     PART 7:  ROTATE THE AIRFOIL BY AN ANGLE ROTANG        *
C        *     ABOUT THE POINT X=ROTCTR.                             *
C        *                                                           *
C        *************************************************************
   35 IF(ROTANG.EQ.0.)  GO TO 37
      SA=SIN(ROTANG)
      CA=COS(ROTANG)
      DO 38 J=JTEBOT,JTETOP
      XX=X(J,1)
      YY=Y(J,1)
      X(J,1)=(XX-ROTCTR)*CA+YY*SA+ROTCTR
   38 Y(J,1)=-(XX-ROTCTR)*SA+YY*CA
C
C        *************************************************************
C        *                                                           *
C        *     PART 8:  MAKE WAKE REGION FOR C-GRID CASES.           *
C        *                                                           *
C        *************************************************************
   37 IF(NTETYP.LT.3.OR.JTEBOT.EQ.1.AND.JTETOP.EQ.JMAX)  GO TO 39
      XT=(X(JTEBOT,1)+X(JTETOP,1))*.5
      YT=(Y(JTEBOT,1)+Y(JTETOP,1))*.5
      XTM=(X(JTEBOT+1,1)+X(JTETOP-1,1))*.5
      DX1=(XT-XTM)*XTFRAC
      DXTOT=XRIGHT-XT
      TALAMR=TAN(ALAMR)
      DELX=.5*TEOPEN*SIN(ALAMR)
      DELY=.5*TEOPEN*COS(ALAMR)
      EL=0.
      DXSUM=0.
      DO 47 J=2,JTEBOT
   47 DXSUM=DXSUM+DX1*(1.+EL*(FLOAT(J-1)**(-WAKEP)))**(J-2)
      FL=DXTOT-DXSUM
      EH=.5
   51 EH=EH*2.
      IF(EH.GE.512.)  GO TO 52
      DXSUM=0.
      DO 44 J=2,JTEBOT
   44 DXSUM=DXSUM+DX1*(1.+EH*(FLOAT(J-1)**(-WAKEP)))**(J-2)
      FH=DXTOT-DXSUM
      IF(FH.GE.0.)  GO TO 51
      FLIM=1.E-7*(XTE-XLE)
C
      DO 40 I=1,1000
      E=EL+FL*(EH-EL)/(FL-FH)
      DXSUM=0.
      DO 45 J=2,JTEBOT
   45 DXSUM=DXSUM+DX1*(1.+E*(FLOAT(J-1)**(-WAKEP)))**(J-2)
      F=DXTOT-DXSUM
      IF(ABS(F).LE.FLIM)  GO TO 48
      IF(F.LT.0..AND.FH.LT.0..OR.F.GT.0..AND.FH.GT.0.)  GO TO 41
      EL=E
      FL=F
      GO TO 40
   41 EH=E
      FH=F
   40 CONTINUE
   52  WRITE(6,102)  EH,I,F
      STOP
C
   48 DX=0.
      DO 42 JJ=2,JTEBOT
      J=JTEBOT+1-JJ
      JJJ=JTETOP-1+JJ
      DX=DX+DX1*(1.+E*(FLOAT(JJ-1)**(-WAKEP)))**(JJ-2)
      X(J,1)=XT+DX-DELX
      X(JJJ,1)=XT+DX+DELX
      Y(J,1)=YT-DX*TALAMR-DELY
   42 Y(JJJ,1)=YT-DX*TALAMR+DELY
C
C        *************************************************************
C        *                                                           *
C        *    PART 9:  IF NDS = 1 CALCULATE DS AS DSI TIMES THE      *
C        *    SURFACE ARCLENGTH.                                     *
C        *                                                           *
C        *************************************************************
   39 IF(NDS.NE.1)  GO TO 46
      JMM=JMAX-1
      JMM2=JMAX-2
      DO 50 J=2,JMM
      ALEN=.5*SQRT((X(J+1,1)-X(J-1,1))**2+(Y(J+1,1)-Y(J-1,1))**2)
   50 DS(J)=ALEN*DSI
      DS(1)=2.*DS(2)-DS(3)
      DS(JMAX)=2.*DS(JMM)-DS(JMM2)
C
   46 IF(JPRT.LT.0.OR.NIBDST.EQ.5.AND.NTETYP.LT.3)  GO TO 2
      GO TO (53,54),NDS
   53 WRITE(6,101)  (J,X(J,1),Y(J,1),DS(J),J=1,JMAX)
      GO TO 55
   54 WRITE(6,100)  (J,X(J,1),Y(J,1),J=1,JMAX)
   55 WRITE(6,110)  TEOPEN
C
    2 RETURN
C
  100 FORMAT(15H0INNER BOUNDARY/3X,1HJ,17X,1HX,17X,1HY/(1X,I3,5X,F13.6,
     1  5X,F13.6))
  101 FORMAT(15H0INNER BOUNDARY/3X,1HJ,17X,1HX,17X,1HY,16X,2HDS/
     1  (1X,I3,5X,F13.6,5X,F13.6,5X,F13.6))
  102 FORMAT(49H0FATAL ERROR IN INNER.  FAILURE TO CONVERGE.  EH=,
     1  F6.1,5X,2HI=,I4,5X,2HF=,F15.12)
  110 FORMAT(8H0TEOPEN=,F13.6)
      END
*DECK INPUT
      SUBROUTINE INPUT
C
*CALL GRIDB
C
*CALL GRIDC
C
*CALL GRID1
C
*CALL GRID2
C
*CALL GRID3
C
      READ(5,GRID1)
      READ(5,GRID2)
      READ(5,GRID3)
      RETURN
      END
*DECK INTERP
      SUBROUTINE INTERP
C
C        THIS SUBROUTINE STARTS WITH THE COARSE GRID SOLUTION
C        AND INTERPOLATES TO FIND THE STARTING SOLUTION FOR THE FINE
C        GRID.
C
*CALL GRIDA
C
*CALL GRIDB
C
*CALL GRIDC
C
      DIMENSION U(170,4),V(170),W(170,2),Z(170)
C
      KMM=KMAX-1
      J2=JS+JINC
      J3=JE-JINC
      J2P=J2+JINC
      J3M=J3-JINC
      JMM=JMAX-1
      JMM2=JMAX-2
C
C     INTERPOLATE X AND Y IN THE ETA(K) DIRECTION.
      K=0
      DO 1 KK=KS,KE,KINC
      K=K+1
    1 V(K)=FLOAT(KK)
C
      DO 2 K=2,KMM
    2 Z(K)=FLOAT(K)
C
      JF=JS
      JL=JE
      IF(NTETYP.LT.3)  GO TO 6
      JF=1
      JL=JMAX
C
    6 DO 3 J=JF,JL,JINC
      K=0
      DO 4 KK=KS,KE,KINC
      K=K+1
      U(K,1)=X(J,KK)
    4 U(K,2)=Y(J,KK)
      CALL CSPLIN(Z,W(1,1),DUM,V,U(1,1),2,KMM,1,K,0)
      CALL CSPLIN(Z,W(1,2),DUM,V,U(1,2),2,KMM,1,K,0)
      DO 5 K=2,KMM
      X(J,K)=W(K,1)
    5 Y(J,K)=W(K,2)
    3 CONTINUE
C
C        INTERPOLATE X AND Y IN THE XI(J) DIRECTION.
      J=1
      DO 11 JJ=JS,JE,JINC
      J=J+1
   11 V(J)=FLOAT(JJ)
      V(1)=FLOAT(JS-JINC)
      J=J+1
      V(J)=FLOAT(JE+JINC)
C
      JF=1
      JL=JMAX
      IF(NTETYP.LT.3)  GO TO 7
      JF=2
      JL=JMM
    7 DO 12 J=JF,JL
   12 Z(J)=FLOAT(J)
C
      DO 13 K=2,KMM
      J=1
      DO 14 JJ=JS,JE,JINC
      J=J+1
      U(J,1)=X(JJ,K)
   14 U(J,2)=Y(JJ,K)
C
      IF(NTETYP.EQ.3)  GO TO 18
      U(1,1)=X(JE,K)
      U(1,2)=Y(JE,K)-TEOPEN
      J=J+1
      U(J,1)=X(JS,K)
      U(J,2)=Y(JS,K)+TEOPEN
      GO TO 19
   18 CONTINUE
      U(1,1)=X(1,K)
      U(1,2)=Y(1,K)
      J=J+1
      U(J,1)=X(JMAX,K)
      U(J,2)=Y(JMAX,K)
   19 CONTINUE
      CALL CSPLIN(Z,X(1,K),DUM,V,U(1,1),JF,JL,1,J,0)
      CALL CSPLIN(Z,Y(1,K),DUM,V,U(1,2),JF,JL,1,J,0)
   13 CONTINUE
C
      IF(NTETYP.GT.1)  GO TO 30
C
C        INTERPOLATE P1, Q1, R1, AND S1 FOR THE CASE OF AN O-GRID
C        WITH A POINT AT THE TRAILING-EDGE.
      J=1
      DO 31 JJ=JS,JE,JINC
      J=J+1
      U(J,1)=P1(JJ)
      U(J,2)=Q1(JJ)
      U(J,3)=R1(JJ)
      U(J,4)=S1(JJ)
   31 V(J)=FLOAT(JJ)
C
      ONE3RD=1./3.
      U(1,1)=P1(3)+(P1(3)-P1(J2))*ONE3RD
      U(1,2)=Q1(3)+(Q1(3)-Q1(J2))*ONE3RD
      U(1,3)=R1(3)+(R1(3)-R1(J2))*ONE3RD
      U(1,4)=S1(3)+(S1(3)-S1(J2))*ONE3RD
      V(1)=1.
      J=J+1
      U(J,1)=P1(JE)+(P1(JE)-P1(J3))*ONE3RD
      U(J,2)=Q1(JE)+(Q1(JE)+Q1(J3))*ONE3RD
      U(J,3)=R1(JE)+(R1(JE)+R1(J3))*ONE3RD
      U(J,4)=S1(JE)+(S1(JE)+S1(J3))*ONE3RD
      V(J)=FLOAT(JMAX)
C
      CALL CSPLIN(Z,P1,DUM,V,U(1,1),1,JMAX,1,J,0)
      CALL CSPLIN(Z,Q1,DUM,V,U(1,2),1,JMAX,1,J,0)
      CALL CSPLIN(Z,R1,DUM,V,U(1,3),1,JMAX,1,J,0)
      CALL CSPLIN(Z,S1,DUM,V,U(1,4),1,JMAX,1,J,0)
      GO TO 21
C
C        INTERPOLATE P1, Q1, R1, AND S1 FOR THE CASE OF AN O-GRID
C        WITH NO POINT AT THE TRAILING EDGE.
   30 IF(NTETYP.EQ.3)  GO TO 20
      DVV=FLOAT((JMAX-5)*JINC)/FLOAT(JE-JS-4*JINC)
      VV=3.-DVV
      J=0
      DO 15 JJ=J2P,J3M,JINC
      J=J+1
      U(J,1)=P1(JJ)
      U(J,2)=Q1(JJ)
      U(J,3)=R1(JJ)
      U(J,4)=S1(JJ)
      VV=VV+DVV
   15 V(J)=VV
      P1(2)=P1(J2)
      P1(JMM)=P1(J3)
      Q1(2)=Q1(J2)
      Q1(JMM)=Q1(J3)
      R1(2)=R1(J2)
      R1(JMM)=R1(J3)
      S1(2)=S1(J2)
      S1(JMM)=S1(J3)
      ONE3RD=1./3.
      R1(1)=R1(J2)+ONE3RD*(R1(J3)-R1(J2))
      S1(1)=S1(J2)+ONE3RD*(S1(J3)-S1(J2))
      TWO3RD=2./3.
      R1(JMAX)=R1(J2)+TWO3RD*(R1(J3)-R1(J2))
      S1(JMAX)=S1(J2)+TWO3RD*(S1(J3)-S1(J2))
C
      CALL CSPLIN(Z,P1,DUM,V,U(1,1),3,JMM2,1,J,0)
      CALL CSPLIN(Z,Q1,DUM,V,U(1,2),3,JMM2,1,J,0)
      CALL CSPLIN(Z,R1,DUM,V,U(1,3),3,JMM2,1,J,0)
      CALL CSPLIN(Z,S1,DUM,V,U(1,4),3,JMM2,1,J,0)
      GO TO 21
C
C        INTERPOLATE P1, Q1, R1, AND S1 FOR THE CASE OF A C-GRID.
   20 J=1
      DO 22 JJ=JS,JE,JINC
      J=J+1
      U(J,1)=P1(JJ)
      U(J,2)=Q1(JJ)
      U(J,3)=R1(JJ)
   22 U(J,4)=S1(JJ)
C
      U(1,1)=-P1(J2)+2.*P1(JS)
      U(1,2)=-Q1(J2)+2.*Q1(JS)
      U(1,3)=-R1(J2)+2.*R1(JS)
      U(1,4)=-S1(J2)+2.*S1(JS)
      J=J+1
      U(J,1)=-P1(J3)+2.*P1(JE)
      U(J,2)=-Q1(J3)+2.*Q1(JE)
      U(J,3)=-R1(J3)+2.*R1(JE)
      U(J,4)=-S1(J3)+2.*S1(JE)
C
      CALL CSPLIN(Z,P1,DUM,V,U(1,1),2,JMM,1,J,0)
      CALL CSPLIN(Z,Q1,DUM,V,U(1,2),2,JMM,1,J,0)
      CALL CSPLIN(Z,R1,DUM,V,U(1,3),2,JMM,1,J,0)
      CALL CSPLIN(Z,S1,DUM,V,U(1,4),2,JMM,1,J,0)
C
   21 RETURN
      END
*DECK OUTER
      SUBROUTINE OUTER
C
C        THIS SUBROUTINE COMPUTES THE LOCATIONS OF POINTS ON THE
C        OUTER BOUNDARY.  IT ALSO FILLS SOME ARRAYS WITH PARAMETERS
C        WHICH DETERMINE THE ANGLES AND SPACING AT THE OUTER BOUNDARY.
C
*CALL GRIDA
C
*CALL GRIDB
C
*CALL GRIDC
C
      DIMENSION PHA(9),SA(13)
C
      DATA N1/1/
C
      YARC(YTARC,X)=YTARC-RARC+SQRT(RARC**2-(X-XARC)**2)
C
C        SPECIAL TREATMENT IF THE USER HAS READ IN HIS OWN OUTER
C        BOUNDARY.
      IF(NOBSHP.LT.6)  GO TO 95
      DO 94 J=1,JMAX
      X(J,KMAX)=XOB(J)
   94 Y(J,KMAX)=YOB(J)
      GO TO 11
C
C        COMPUTE LOCATIONS OF POINTS ON OUTER BOUNDARY FOR THE CASE
C        OF AN O-GRID AND THE OUTER BOUNDARY BEING A CIRCLE (NOBSHP=1).
   95 HTEOPN=.5*TEOPEN
      IF(NOBSHP.GT.1)  GO TO 2
      PHTEO=ASIN(HTEOPN/RADOB)
      DTH=(TWOPI-2.*PHTEO)/FLOAT(JMAX)
      TH=-DTH+PHTEO
      IF(NTETYP.GT.1)  TH=-DTH*.5+PHTEO
      TH=TH+ALAMR
      IF(JPRT.GE.0)  WRITE(6,101)
C
      DO 10 J=1,JMAX
      TH=TH+DTH
      PH=TH
      IF(NOBDST.EQ.3)  PH=OBANGS(J)
      X(J,KMAX)=RADOB*COS(PH)+XOBCNT
      Y(J,KMAX)=-RADOB*SIN(PH)
      IF(JPRT.LT.0)  GO TO 10
      PHD=PH*TODEG
      WRITE(6,100)  J,PHD,X(J,KMAX),Y(J,KMAX)
   10 CONTINUE
      GO TO 11
C
C        PREPARE FOR COMPUTING RECTANGULAR OUTER BOUNDARIES.
    2 IF(NOBSHP.EQ.3.OR.NOBSHP.EQ.5)  GO TO 60
      YTMRC=YTOP-RCORN
      YBPRC=YBOTOM+RCORN
      XRMRC=XRIGHT-RCORN
      XLPRC=XLEFT+RCORN
      RRCORN=1.0/RCORN
      QCIRC=.5*PI*RCORN
      IF(NOBSHP.EQ.4)  GO TO 4
      IF(NOBDST.GT.1)  GO TO 12
C
C        COMPUTE LOCATIONS OF POINTS ON OUTER BOUNDARY FOR THE CASE
C        OF AN O-GRID, AND THE OUTER BOUNDARY BEING A RECTANGLE WITH
C        ROUNDED CORNERS, AND THE POINTS DISTRIBUTED BY EQUAL ARC-
C        LENGTH INCREMENTS (NOBSHP=2, NOBDST=1).
      STOTAL=2.*(XRIGHT-XLEFT)+2.*(YTOP-YBOTOM)-8.*RCORN+4.*QCIRC-TEOPEN
      DA=STOTAL/FLOAT(JMAX)
      SA(1)=-YBPRC-HTEOPN
      SA(2)=SA(1)+QCIRC
      SA(3)=SA(2)+XRMRC-XLPRC
      SA(4)=SA(3)+QCIRC
      SA(5)=SA(4)+YTMRC-YBPRC
      SA(6)=SA(5)+QCIRC
      SA(7)=SA(6)+XRMRC-XLPRC
      SA(8)=SA(7)+QCIRC
      SA(9)=STOTAL
      S=0.
      X(1,KMAX)=XRIGHT
      Y(1,KMAX)=0.-HTEOPN
      IF(NTETYP.EQ.0)  GO TO 16
      S=.5*DA
      Y(1,KMAX)=-S-HTEOPN
   16 KEY=1
      IF(JPRT.GE.0)  WRITE(6,105)
      IF(JPRT.GE.0)  WRITE(6,104) N1,KEY,S,X(1,KMAX),Y(1,KMAX)
C
      DO 15 J=2,JMAX
      S=S+DA
    9 IF(S.LE.SA(KEY))  GO TO 8
      KEY=KEY+1
      GO TO 9
C
    8 IF(KEY.GT.1)  CS=S-SA(KEY-1)
      GO TO (31,32,33,34,35,36,37,38,39),KEY
   31 X(J,KMAX)=XRIGHT
      Y(J,KMAX)=Y(J-1,KMAX)-DA
      GO TO 17
   32 PHIB=-CS*RRCORN
      X(J,KMAX)=XRMRC+RCORN*COS(PHIB)
      Y(J,KMAX)=YBPRC+RCORN*SIN(PHIB)
      GO TO 17
   33 X(J,KMAX)=XRMRC-CS
      Y(J,KMAX)=YBOTOM
      GO TO 17
   34 PHIB=-PIO2-CS*RRCORN
      X(J,KMAX)=XLPRC+RCORN*COS(PHIB)
      Y(J,KMAX)=YBPRC+RCORN*SIN(PHIB)
      GO TO 17
   35 X(J,KMAX)=XLEFT
      Y(J,KMAX)=YBPRC+CS
      GO TO 17
   36 PHIB=PI-CS*RRCORN
      X(J,KMAX)=XLPRC+RCORN*COS(PHIB)
      Y(J,KMAX)=YTMRC+RCORN*SIN(PHIB)
      GO TO 17
   37 X(J,KMAX)=XLPRC+CS
      Y(J,KMAX)=YTOP
      GO TO 17
   38 PHIB=PIO2-CS*RRCORN
      X(J,KMAX)=XRMRC+RCORN*COS(PHIB)
      Y(J,KMAX)=YTMRC+RCORN*SIN(PHIB)
      GO TO 17
   39 X(J,KMAX)=XRIGHT
      Y(J,KMAX)=YTMRC-CS
   17 CONTINUE
      IF(JPRT.LT.0)  GO TO 15
      WRITE(6,104)  J,KEY,S,X(J,KMAX),Y(J,KMAX)
   15 CONTINUE
      GO TO 11
C
C        COMPUTE LOCATIONS OF POINTS ON OUTER BOUNDARY FOR THE CASE
C        OF AN O-GRID, AND THE OUTER BOUNDARY BEING A RECTANGLE WITH
C        ROUNDED CORNERS, AND THE POINTS DISTRIBUTED IN SOME ANGULAR
C        FASHION (NOBSHP=2, NOBDST=2 OR 3).
   12 PHTEO=ATAN2(HTEOPN,XRIGHT)
      DTH=TWOPI/FLOAT(JMAX)
      TH=-DTH+PHTEO
      IF(NTETYP.GT.1)  TH=-DTH*.5
      KEY=1
      PHA(1)=-ATAN2(YBPRC,XRIGHT)
      PHA(2)=-ATAN2(YBOTOM,XRMRC)
      PHA(3)=-ATAN2(YBOTOM,XLPRC)
      PHA(4)=-ATAN2(YBPRC,XLEFT)
      PHA(5)=TWOPI-ATAN2(YTMRC,XLEFT)
      PHA(6)=TWOPI-ATAN2(YTOP,XLPRC)
      PHA(7)=TWOPI-ATAN2(YTOP,XRMRC)
      PHA(8)=TWOPI-ATAN2(YTMRC,XRIGHT)
      PHA(9)=TWOPI
      IF(JPRT.GE.0)  WRITE(6,102)
C
      DO 13 J=1,JMAX
      TH=TH+DTH
      PH=TH
      IF(NOBDST.EQ.3)  PH=OBANGS(J)
   20 IF(PH.LE.PHA(KEY).OR.KEY.EQ.9)  GO TO 19
      KEY=KEY+1
      GO TO 20
C
   19 GO TO (21,22,23,24,25,26,27,28,29),KEY
   21 X(J,KMAX)=XRIGHT
      Y(J,KMAX)=-XRIGHT*TAN(PH)
      GO TO 14
   22 PHIB=-PH+ASIN((-XRMRC*TAN(PH)-YBPRC)*COS(PH)*RRCORN)
      X(J,KMAX)=XRMRC+RCORN*COS(PHIB)
      Y(J,KMAX)=YBPRC+RCORN*SIN(PHIB)
      GO TO 14
   23 X(J,KMAX)=-YBOTOM/TAN(PH)
      Y(J,KMAX)=YBOTOM
      GO TO 14
   24 PHIB=-PH+ASIN((-XLPRC*TAN(PH)-YBPRC)*COS(PH)*RRCORN)
      X(J,KMAX)=XLPRC+RCORN*COS(PHIB)
      Y(J,KMAX)=YBPRC+RCORN*SIN(PHIB)
      GO TO 14
   25 X(J,KMAX)=XLEFT
      Y(J,KMAX)=-XLEFT*TAN(PH)
      GO TO 14
   26 PHIB=-PH+ASIN((-XLPRC*TAN(PH)-YTMRC)*COS(PH)*RRCORN)
      X(J,KMAX)=XLPRC+RCORN*COS(PHIB)
      Y(J,KMAX)=YTMRC+RCORN*SIN(PHIB)
      GO TO 14
   27 X(J,KMAX)=-YTOP/TAN(PH)
      Y(J,KMAX)=YTOP
      GO TO 14
   28 PHIB=-PH+ASIN((-XRMRC*TAN(PH)-YTMRC)*COS(PH)*RRCORN)
      X(J,KMAX)=XRMRC+RCORN*COS(PHIB)
      Y(J,KMAX)=YTMRC+RCORN*SIN(PHIB)
      GO TO 14
   29 X(J,KMAX)=XRIGHT
      Y(J,KMAX)=-XRIGHT*TAN(PH)
   14 CONTINUE
      IF(JPRT.LT.0)  GO TO 13
      PHD=PH*TODEG
      WRITE(6,103)  J,KEY,PHD,X(J,KMAX),Y(J,KMAX)
   13 CONTINUE
      GO TO 11
C
C        PREPARE FOR COMPUTING OUTER BOUNDARY FOR CASCADE GRIDS.
   60 XRMRC=XRIGHT-RCORN
      XLPRC=XLEFT+RCORN
      RRCORN=1./RCORN
      SALR=SIN(ALAMR)
      CALR=COS(ALAMR)
      TALR=TAN(ALAMR)
      SALF=SIN(ALAMF)
      CALF=COS(ALAMF)
      TALF=TAN(ALAMF)
      XARC=(XLE*SALR+XTE*SALF)/(SALR+SALF)
      IF(SALR.NE.0.)  RARC=(XTE-XARC)/SALR
      IF(SALF.NE.0.)  RARC=(XARC-XLE)/SALF
      RRARC=1./RARC
C
      Y10=YARC(YTOP,XTE)
      X11=XRMRC+RCORN*COS(PIO2-ALAMR)
      Y11=Y10-(X11-XTE)*TALR
      Y12=Y11-RCORN*SIN(PIO2-ALAMR)
      Y13=-(XRIGHT-XTE)*TALR
      Y14=Y10-(XRIGHT-XTE)*TALR
C
      Y3=YARC(YBOTOM,XTE)
      X2=XRIGHT-RCORN*(SALR+1.)
      Y2=Y3-(X2-XTE)*TALR
      Y1=Y2+RCORN*SIN(PIO2-ALAMR)
      Y15=Y3-(XRIGHT-XTE)*TALR
      Y16=Y11+(X11-X2)*TALR
C
      Y4=YARC(YBOTOM,XLE)
      X5=XLEFT+RCORN*(SALF+1.)
      Y5=Y4-(XLE-X5)*TALF
      Y6=Y5+RCORN*SIN(PIO2-ALAMF)
C
      Y9=YARC(YTOP,XLE)
      X8=XLPRC-RCORN*COS(PIO2-ALAMF)
      Y8=Y9-(XLE-X8)*TALF
      Y7=Y8-RCORN*SIN(PIO2-ALAMF)
C
      IF(JPRT.GE.0)  WRITE(6,106)  XRIGHT,Y1,X2,Y2,XTE,Y3,XLE,Y4,
     1  X5,Y5,XLEFT,Y6,XLEFT,Y7,X8,Y8,XLE,Y9,XTE,Y10,X11,Y11,XRIGHT,
     2  Y12,XRIGHT,Y13,XRIGHT,Y14,XRIGHT,Y15,X2,Y16
      IF(NOBSHP.EQ.5)  GO TO 80
C
C        COMPUTE LOCATIONS OF OUTER BOUNDARY POINTS FOR THE CASE OF AN
C        O-GRID FOR A CASCADE, AND AN EQUAL-ARC-LENGTH DISTRIBUTION
C        (NOBSHP=3).
      DS3T=SQRT((X2-XTE)**2+(Y2-Y3)**2)
      DS5T=SQRT((XLE-X5)**2+(Y4-Y5)**2)
      DS9T=SQRT((XLE-X8)**2+(Y9-Y8)**2)
      DS16T=SQRT((X2-XTE)**2+(Y16-Y10)**2)
      DS11T=SQRT((X11-X2)**2+(Y11-Y16)**2)
      STOT1=Y13-Y1+(PIO2+ALAMR)*RCORN-HTEOPN
      STOT2=Y7-Y6+DS3T+DS5T+DS9T+DS16T+PI*RCORN+2.*RARC*(ALAMF+ALAMR)
      STOT3=DS11T+(PIO2-ALAMR)*RCORN+Y12-Y13-HTEOPN
      STOTAL=STOT1+STOT2+STOT3
      RST=1./STOTAL
      FJ1=FLOAT(JMAX)*STOT1*RST
      FJ3=FLOAT(JMAX)*STOT3*RST
      IF(NTETYP.EQ.1)  GO TO 96
      FJ1=FJ1-0.5
      FJ3=FJ3+0.5
   96 JN1=IFIX(FJ1+0.5)
      JN3=IFIX(FJ3+0.5)
      JN2=JMAX-JN1-JN3
      FJN1=FLOAT(JN1)
      FJN3=FLOAT(JN3)
      IF(NTETYP.EQ.1)  GO TO 97
      FJN1=FJN1+0.5
      FJN3=FJN3-0.5
   97 DA1=STOT1/FJN1
      DA2=STOT2/FLOAT(JN2)
      DA3=STOT3/FJN3
      JN1P=JN1+1
      JN12P=JN1+JN2+1
      SA(1)=Y13-Y1-HTEOPN
      SA(2)=SA(1)+RCORN*(PIO2+ALAMR)
      SA(3)=SA(2)+DS3T
      SA(4)=SA(3)+RARC*(ALAMF+ALAMR)
      SA(5)=SA(4)+DS5T
      SA(6)=SA(5)+RCORN*(PIO2+ALAMF)
      SA(7)=SA(6)+Y7-Y6
      SA(8)=SA(7)+RCORN*(PIO2-ALAMF)
      SA(9)=SA(8)+DS9T
      SA(10)=SA(9)+RARC*(ALAMF+ALAMR)
      SA(11)=SA(10)+DS11T+DS16T
      SA(12)=SA(11)+RCORN*(PIO2-ALAMR)
      SA(13)=STOTAL
      RDS2=1./(SA(2)-SA(1))
      RDS4=1./(SA(4)-SA(3))
      RDS6=1./(SA(6)-SA(5))
      RDS8=1./(SA(8)-SA(7))
      RDS10=1./(SA(10)-SA(9))
      RDS12=1./(SA(12)-SA(11))
      S=0.
      X(1,KMAX)=XRIGHT
      Y(1,KMAX)=Y13-HTEOPN
      IF(NTETYP.EQ.1)  GO TO 3
      S=0.5*DA1
      Y(1,KMAX)=Y13-S-HTEOPN
    3 KEY=1
      IF(JPRT.GE.0) WRITE(6,107)
      THETOB(1)=PIO2
      THD=THETOB(1)*TODEG
      IF(JPRT.GE.0)  WRITE(6,104)  N1,KEY,S,X(1,KMAX),Y(1,KMAX),THD
C
      DO 74 J=2,JMAX
      IF(J.LE.JN1P)  S=S+DA1
      IF(J.GT.JN1P.AND.J.LE.JN12P)  S=S+DA2
      IF(J.GT.JN12P)  S=S+DA3
   79 IF(S.LE.SA(KEY)) GO TO 78
      KEY=KEY+1
      GO TO 79
   78 IF(KEY.GT.1)  CS=S-SA(KEY-1)
C
      GO TO (61,62,63,64,65,66,67,68,69,70,71,72,73),KEY
   61 X(J,KMAX)=XRIGHT
      Y(J,KMAX)=Y(J-1,KMAX)-DA1
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2
      GO TO 75
   62 PHIB=-CS*RRCORN
      X(J,KMAX)=XRMRC+RCORN*COS(PHIB)
      Y(J,KMAX)=Y1+RCORN*SIN(PHIB)
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2-ALAMR*CS*RDS2
      GO TO 75
   63 X(J,KMAX)=X2-CS*CALR
      Y(J,KMAX)=Y2+CS*SALR
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2-ALAMR
      GO TO 75
   64 PHIB=PIO2-ALAMR+CS*RRARC
      X(J,KMAX)=XARC+RARC*COS(PHIB)
      Y(J,KMAX)=YARC(YBOTOM,X(J,KMAX))
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2-ALAMR+(ALAMF+ALAMR)*CS*RDS4
      GO TO 75
   65 X(J,KMAX)=XLE-CS*CALF
      Y(J,KMAX)=Y4-CS*SALF
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2+ALAMF
      JCORN(2)=J
      GO TO 75
   66 PHIB=-PIO2+ALAMF-CS*RRCORN
      X(J,KMAX)=XLPRC+RCORN*COS(PHIB)
      Y(J,KMAX)=Y6+RCORN*SIN(PHIB)
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2+ALAMF*(1.-CS*RDS6)
      GO TO 75
   67 X(J,KMAX)=XLEFT
      Y(J,KMAX)=Y6+CS
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2
      GO TO 75
   68 PHIB=PI-CS*RRCORN
      X(J,KMAX)=XLPRC+RCORN*COS(PHIB)
      Y(J,KMAX)=Y7+RCORN*SIN(PHIB)
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2+ALAMF*CS*RDS8
      GO TO 75
   69 X(J,KMAX)=X8+CS*CALF
      Y(J,KMAX)=Y8+CS*SALF
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2+ALAMF
      GO TO 75
   70 PHIB=PIO2+ALAMF-CS*RRARC
      X(J,KMAX)=XARC+RARC*COS(PHIB)
      Y(J,KMAX)=YARC(YTOP,X(J,KMAX))
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2+ALAMF-(ALAMR+ALAMF)*CS*RDS10
      GO TO 75
   71 X(J,KMAX)=XTE+CS*CALR
      Y(J,KMAX)=Y10-CS*SALR
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2-ALAMR
      GO TO 75
   72 PHIB=PIO2-CS*RRCORN-ALAMR
      X(J,KMAX)=XRMRC+RCORN*COS(PHIB)
      Y(J,KMAX)=Y12+RCORN*SIN(PHIB)
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2-ALAMR*(1.-CS*RDS12)
      GO TO 75
   73 X(J,KMAX)=XRIGHT
      Y(J,KMAX)=Y12-CS
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2
   75 CONTINUE
      IF(JPRT.LT.0)  GO TO 74
      THD=THETOB(J)*TODEG
      WRITE(6,104)  J,KEY,S,X(J,KMAX),Y(J,KMAX),THD
   74 CONTINUE
      JCORN(3)=2*JN1P+JN2-JCORN(2)
      JCORN(1)=JN1P
      JCORN(4)=JN1P+JN2
      WRITE(6,108)  JCORN
      GO TO 11
C
C        COMPUTE LOCATIONS OF POINTS ON OUTER BOUNDARY FOR THE CASE
C        OF A C-GRID, AND THE OUTER BOUNDARY BEING A RECTANGLE WITH
C        ROUNDED CORNERS ON THE LEFT, AND THE POINTS DISTRIBUTED
C        BY EQUAL ARC-LENGTH INCREMENTS (NOBSHP=4, NOBDST=1).
    4 IF(NOBDST.GT.1)  GO TO 1
      STOTAL=2.*(XRIGHT-XLEFT)+YTOP-YBOTOM-4.*RCORN+2.*QCIRC
      DA=STOTAL/FLOAT(JMAX-1)
      SA(1)=XRIGHT-XLPRC
      SA(2)=SA(1)+QCIRC
      SA(3)=SA(2)+YTMRC-YBPRC
      SA(4)=SA(3)+QCIRC
      SA(5)=SA(4)+XRIGHT-XLPRC
      X(1,KMAX)=XRIGHT
      Y(1,KMAX)=YBOTOM
      S=0.
      KEY=1
      IF(JPRT.GE.0)  WRITE(6,105)
      IF(JPRT.GE.0)  WRITE(6,104) N1,KEY,S,X(1,KMAX),Y(1,KMAX)
C
      DO 40 J=2,JMAX
      S=S+DA
    6 IF(S.LE.SA(KEY))  GO TO 7
      KEY=KEY+1
      GO TO 6
C
    7 IF(KEY.GT.1)  CS=S-SA(KEY-1)
      GO TO (41,42,43,44,45),KEY
   41 X(J,KMAX)=X(J-1,KMAX)-DA
      Y(J,KMAX)=YBOTOM
      GO TO 46
   42 PHIB=-PIO2-CS*RRCORN
      X(J,KMAX)=XLPRC+RCORN*COS(PHIB)
      Y(J,KMAX)=YBPRC+RCORN*SIN(PHIB)
      GO TO 46
   43 X(J,KMAX)=XLEFT
      Y(J,KMAX)=YBPRC+CS
      GO TO 46
   44 PHIB=PI-CS*RRCORN
      X(J,KMAX)=XLPRC+RCORN*COS(PHIB)
      Y(J,KMAX)=YTMRC+RCORN*SIN(PHIB)
      GO TO 46
   45 X(J,KMAX)=XLPRC+CS
      Y(J,KMAX)=YTOP
   46 CONTINUE
      IF(JPRT.LT.0)  GO TO 40
      WRITE(6,104)  J,KEY,S,X(J,KMAX),Y(J,KMAX)
   40 CONTINUE
      GO TO 11
C
C        COMPUTE LOCATIONS OF POINTS ON OUTER BOUNDARY FOR THE CASE
C        OF A C-GRID, AND THE OUTER BOUNDARY BEING A RECTANGLE WITH
C        ROUNDED CORNERS ON THE LEFT, AND THE POINTS DISTRIBUTED
C        IN SOME ANGULAR FASHION (NOBSHP=4, NOBDST=2 OR 3).
    1 THT=-ATAN2(YTOP,XRIGHT)
      THB=-ATAN2(YBOTOM,XRIGHT)
      THTOT=TWOPI+THT-THB
      DTH=THTOT/FLOAT(JMAX-1)
      TH=THB-DTH
      KEY=1
      PHA(1)=-ATAN2(YBOTOM,XLPRC)
      PHA(2)=-ATAN2(YBPRC,XLEFT)
      PHA(3)=TWOPI-ATAN2(YTMRC,XLEFT)
      PHA(4)=TWOPI-ATAN2(YTOP,XLPRC)
      PHA(5)=TWOPI+THT
      IF(JPRT.GE.0.AND.NOBSHP.EQ.4)  WRITE(6,102)
C
      DO 50 J=1,JMAX
      TH=TH+DTH
      PH=TH
      IF(NOBDST.EQ.3)  PH=OBANGS(J)
   56 IF(PH.LE.PHA(KEY).OR.KEY.EQ.5)  GO TO 57
      KEY=KEY+1
      GO TO 56
C
   57 GO TO (51,52,53,54,55),KEY
   51 X(J,KMAX)=-YBOTOM/TAN(PH)
      Y(J,KMAX)=YBOTOM
      GO TO 58
   52 PHIB=-PH+ASIN((-XLPRC*TAN(PH)-YBPRC)*COS(PH)*RRCORN)
      X(J,KMAX)=XLPRC+RCORN*COS(PHIB)
      Y(J,KMAX)=YBPRC+RCORN*SIN(PHIB)
      GO TO 58
   53 X(J,KMAX)=XLEFT
      Y(J,KMAX)=-XLEFT*TAN(PH)
      GO TO 58
   54 PHIB=-PH+ASIN((-XLPRC*TAN(PH)-YTMRC)*COS(PH)*RRCORN)
      X(J,KMAX)=XLPRC+RCORN*COS(PHIB)
      Y(J,KMAX)=YTMRC+RCORN*SIN(PHIB)
      GO TO 58
   55 X(J,KMAX)=-YTOP/TAN(PH)
      Y(J,KMAX)=YTOP
   58 CONTINUE
      IF(JPRT.LT.0)  GO TO 50
      PHD=PH*TODEG
      WRITE(6,103)  J,KEY,PHD,X(J,KMAX),Y(J,KMAX)
   50 CONTINUE
      GO TO 11
C
C        COMPUTE LOCATIONS OF OUTER BOUNDARY POINTS FOR THE CASE OF A
C        C-GRID FOR A CASCADE (NOBSHP=5).
   80 DS1T=SQRT((XRIGHT-XTE)**2+(Y15-Y3)**2)
      DS3T=SQRT((XLE-X5)**2+(Y4-Y5)**2)
      DS7T=SQRT((XLE-X8)**2+(Y9-Y8)**2)
      DS9T=SQRT((XRIGHT-XTE)**2+(Y14-Y10)**2)
      STOTAL=DS1T+DS3T+DS7T+DS9T+PI*RCORN+2.*RARC*(ALAMF+ALAMR)+Y7-Y6
      DA=STOTAL/FLOAT(JMAX-1)
      SA(1)=DS1T
      SA(2)=SA(1)+RARC*(ALAMF+ALAMR)
      SA(3)=SA(2)+DS3T
      SA(4)=SA(3)+RCORN*(PIO2+ALAMF)
      SA(5)=SA(4)+Y7-Y6
      SA(6)=SA(5)+RCORN*(PIO2-ALAMF)
      SA(7)=SA(6)+DS7T
      SA(8)=SA(7)+RARC*(ALAMF+ALAMR)
      SA(9)=STOTAL
      RDS2=1./(SA(2)-SA(1))
      RDS4=1./(SA(4)-SA(3))
      RDS6=1./(SA(6)-SA(5))
      RDS8=1./(SA(8)-SA(7))
      X(1,KMAX)=XRIGHT
      Y(1,KMAX)=Y15
      S=0.
      KEY=1
      IF(JPRT.GE.0)  WRITE(6,107)
      THETOB(1)=PIO2-ALAMR
      THD=THETOB(1)*TODEG
      IF(JPRT.GE.0)  WRITE(6,104)  N1,KEY,S,X(1,KMAX),Y(1,KMAX),THD
C
      DO 90 J=2,JMAX
      S=S+DA
   91 IF(S.LE.SA(KEY))  GO TO 92
      KEY=KEY+1
      GO TO 91
   92 IF(KEY.GT.1)  CS=S-SA(KEY-1)
C
      GO TO (81,82,83,84,85,86,87,88,89),KEY
   81 X(J,KMAX)=XRIGHT-S*CALR
      Y(J,KMAX)=Y15+S*SALR
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2-ALAMR
      GO TO 93
   82 PHIB=PIO2-ALAMR+CS*RRARC
      X(J,KMAX)=XARC+RARC*COS(PHIB)
      Y(J,KMAX)=YARC(YBOTOM,X(J,KMAX))
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2-ALAMR+(ALAMF+ALAMR)*CS*RDS2
      GO TO 93
   83 X(J,KMAX)=XLE-CS*CALF
      Y(J,KMAX)=Y4-CS*SALF
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2+ALAMF
      JCORN(2)=J
      GO TO 93
   84 PHIB=-PIO2+ALAMF-CS*RRCORN
      X(J,KMAX)=XLPRC+RCORN*COS(PHIB)
      Y(J,KMAX)=Y6+RCORN*SIN(PHIB)
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2+ALAMF*(1.-CS*RDS4)
      GO TO 93
   85 X(J,KMAX)=XLEFT
      Y(J,KMAX)=Y6+CS
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2
      GO TO 93
   86 PHIB=PI-CS*RRCORN
      X(J,KMAX)=XLPRC+RCORN*COS(PHIB)
      Y(J,KMAX)=Y7+RCORN*SIN(PHIB)
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2+ALAMF*CS*RDS6
      GO TO 93
   87 X(J,KMAX)=X8+CS*CALF
      Y(J,KMAX)=Y8+CS*SALF
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2+ALAMF
      GO TO 93
   88 PHIB=PIO2+ALAMF-CS*RRARC
      X(J,KMAX)=XARC+RARC*COS(PHIB)
      Y(J,KMAX)=YARC(YTOP,X(J,KMAX))
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2+ALAMF-(ALAMR+ALAMF)*CS*RDS8
      GO TO 93
   89 X(J,KMAX)=XTE+CS*CALR
      Y(J,KMAX)=Y10-CS*SALR
      IF(NOBCAS.EQ.0)  THETOB(J)=PIO2-ALAMR
   93 CONTINUE
      IF(JPRT.LT.0)  GO TO 90
      THD=THETOB(J)*TODEG
      WRITE(6,104)  J,KEY,S,X(J,KMAX),Y(J,KMAX),THD
   90 CONTINUE
      JCORN(1)=1
      JCORN(4)=JMAX
      JCORN(3)=JMAX+1-JCORN(2)
      WRITE(6,108) JCORN
      GO TO 11
C
   11 CONTINUE
      RETURN
C
  100 FORMAT(1X,I3,5X,F8.2,2(5X,F13.6))
  101 FORMAT(15H0OUTER BOUNDARY/3X,1HJ,10X,3HPHI,17X,1HX,17X,1HY)
  102 FORMAT(15H0OUTER BOUNDARY/3X,1HJ,4X,3HKEY,10X,3HPHI,17X,1HX,17X,
     1  1HY)
  103 FORMAT(1X,I3,5X,I2,5X,F8.2,2(5X,F13.6))
  104 FORMAT(1X,I3,5X,I2,4(5X,F13.6))
  105 FORMAT(15H0OUTER BOUNDARY/3X,1HJ,4X,3HKEY,9X,9HARCLENGTH,17X,1HX,
     1  17X,1HY)
  106 FORMAT(49H0OUTER BOUNDARY FOR CASCADE GRIDS IS COMPOSED OF ,
     1  39HSTRAIGHT LINES JOINED BY CIRCULAR ARCS./19H TRANSITION POINTS
     * ,
     2  29HBETWEEN LINE SEGMENTS FOLLOW./13X,1HX,17X,1HY/(1X,F13.6,
     3  5X,F13.6))
  107 FORMAT(1H0,2X,1HJ,4X,3HKEY,9X,9HARCLENGTH,17X,1HX,17X,1HY,12X,
     1  6HTHETOB)
  108 FORMAT(7H0JCORN=,4I10)
      END
*DECK OUTPUT
      SUBROUTINE OUTPUT
C
C        THIS SUBROUTINE DOES ALL OF THE OUTPUT OF THE FINISHED
C        GRID, BY PRINTING ON UNIT 6, AND BY UNFORMATTED WRITES
C        ON UNIT 7.
C
*CALL GRIDA
C
*CALL GRIDB
C
*CALL GRIDC
C
      COMMON/SCRACH/AJAC(85),XYMETS(85,4),DUM(1021)
C
      IF(JPRT.LE.0.AND.NOUT.EQ.0)  GO TO 4
C
      IF(JPRT.GT.0)  WRITE(6,112)
      IF(NOUT.GT.0.AND.NTETYP.LT.3)  WRITE(7)  JMAX,KMAX,TEOPEN,JCORN
      IF(NOUT.GT.0.AND.NTETYP.EQ.3)  WRITE(7)  JTEBOT,JTETOP,JMAX,KMAX,
     1  TEOPEN,JCORN
C
C        FOR CASCADE CASES, THE UPPER AND LOWER OUTER BOUNDARY POINTS
C        LINE UP FOR J BETWEEN JCORN(1) AND JCORN(2), INCLUSIVE,
C        ON THE BOTTOM, AND FOR J BETWEEN JCORN(3) AND JCORN (4),
C        INCLUSIVE, ON THE TOP.
C
      DO 26 J=1,JMAX
      IF(NOUT.LT.2)  GO TO 1
C
      DO 24 K=1,KMAX
      IF(J.NE.1)  GO TO 20
      XX=0.5*(-3.*X(J,K)+4.*X(J+1,K)-X(J+2,K))
      YX=0.5*(-3.*Y(J,K)+4.*Y(J+1,K)-Y(J+2,K))
      GO TO 21
   20 IF(J.EQ.JMAX)  GO TO 22
      XX=0.5*(-X(J-1,K)+X(J+1,K))
      YX=0.5*(-Y(J-1,K)+Y(J+1,K))
      GO TO 21
   22 XX=0.5*(X(J-2,K)-4.*X(J-1,K)+3.*X(J,K))
      YX=0.5*(Y(J-2,K)-4.*Y(J-1,K)+3.*Y(J,K))
   21 IF(K.NE.1)  GO TO 23
      XE=0.5*(-3.*X(J,K)+4.*X(J,K+1)-X(J,K+2))
      YE=0.5*(-3.*Y(J,K)+4.*Y(J,K+1)-Y(J,K+2))
      GO TO 2
   23 IF(K.EQ.KMAX)  GO TO 25
      XE=0.5*(-X(J,K-1)+X(J,K+1))
      YE=0.5*(-Y(J,K-1)+Y(J,K+1))
      GO TO 2
   25 XE=0.5*(X(J,K-2)-4.*X(J,K-1)+3.*X(J,K))
      YE=0.5*(Y(J,K-2)-4.*Y(J,K-1)+3.*Y(J,K))
C
    2 AJAC(K)=XX*YE-XE*YX
      IF(NOUT.EQ.2)  GO TO 24
C
      XYMETS(K,1)=XX
      XYMETS(K,2)=XE
      XYMETS(K,3)=YX
      XYMETS(K,4)=YE
   24 CONTINUE
C
    1 IF(NOUT.EQ.0)  GO TO 3
      GO TO (11,12,13),NOUT
   11 WRITE(7)  (X(J,K),K=1,KMAX),(Y(J,K),K=1,KMAX)
      IF(J.EQ.JMAX)  WRITE(6,113)
      GO TO 3
   12 WRITE(7)  (X(J,K),K=1,KMAX),(Y(J,K),K=1,KMAX),(AJAC(K),K=1,KMAX)
      IF(J.EQ.JMAX)  WRITE(6,114)
      GO TO 3
   13 WRITE(7)  (X(J,K),K=1,KMAX),(Y(J,K),K=1,KMAX),(AJAC(K),K=1,KMAX),
     1  ((XYMETS(K,L),K=1,KMAX),L=1,4)
      IF(J.EQ.JMAX)  WRITE(6,115)
C
    3 IF(JPRT.LE.0)  GO TO 26
      IF(((J-1)/JPRT)*JPRT.NE.J-1)  GO TO 26
C
      WRITE(6,100)  J
      WRITE(6,116)  (X(J,K),K=1,KMAX)
      WRITE(6,101)  J
      WRITE(6,116)  (Y(J,K),K=1,KMAX)
C
      IF(NOUT.LT.2)  GO TO 26
      WRITE(6,102)  J
      WRITE(6,116)  (AJAC(K),K=1,KMAX)
C
      IF(NOUT.EQ.2)  GO TO 26
      WRITE(6,103)  J
      WRITE(6,116)  (XYMETS(K,1),K=1,KMAX)
      WRITE(6,104)  J
      WRITE(6,116)  (XYMETS(K,2),K=1,KMAX)
      WRITE(6,105)  J
      WRITE(6,116)  (XYMETS(K,3),K=1,KMAX)
      WRITE(6,106)  J
      WRITE(6,116)  (XYMETS(K,4),K=1,KMAX)
   26 CONTINUE
C
    4 RETURN
C
  100 FORMAT(8H0X AT J=,I3)
  101 FORMAT(8H0Y AT J=,I3)
  102 FORMAT(16H0JACOBIANS AT J=,I3)
  103 FORMAT(13H0DX/DXI AT J=,I3)
  104 FORMAT(14H0DX/DETA AT J=,I3)
  105 FORMAT(13H0DY/DXI AT J=,I3)
  106 FORMAT(14H0DY/DETA AT J=,I3)
  116 FORMAT(1X,8E15.6)
  112 FORMAT(1H1,50X,13HFINISHED GRID/)
  113 FORMAT(50H0SOLUTION VALUES FOR X AND Y HAVE BEEN WRITTEN ON ,
     1  7HUNIT 7.)
  114 FORMAT(40H0SOLUTION VALUES FOR X, Y, AND JACOBIAN ,
     1  28HHAVE BEEN WRITTEN ON UNIT 7.)
  115 FORMAT(50H0SOLUTION VALUES FOR X, Y, JACOBIAN, AND THE FOUR ,
     1  46HMETRIC QUANTITIES HAVE BEEN WRITTEN ON UNIT 7.)
      END
*DECK RELAX
      SUBROUTINE RELAX
C
C        THIS SUBROUTINE APPLIES AN S.L.O.R. TECHNIQUE TO SOLVE
C        THE POISSON EQUATIONS THAT GENERATE THE GRID.
C
*CALL GRIDA
C
*CALL GRIDB
C
*CALL GRIDC
C
      COMMON/SCRACH/RHS1(170),RHS2(170),RHS3(170),RHS4(170),PCORCT(170),
     1  PRESID(170),QRESID(170),RRESID(170),DUMMY(86)
C
      DIMENSION A(170),B(170),C(170),D(170),F(170),H(170),JPPER(170),
     1  JMPER(170),XXB(170),YXB(170),XXXB(170),F0(170),F1(170),F2(170),
     2  YXXB(170),XEB(170),YEB(170),XXC(170),YXC(170),XXXC(170),
     3  YXXC(170),XEC(170),YEC(170),F3(170),F4(170),F5(170),RJC(170),
     4  RJB(170),G(170),S(170),QPS(170),XF(170),YG(170),XC(170),
     5  RCORCT(170),SCORCT(170),XXD(170),YXD(170),XED(170),
     6  YED(170),AJSQ(170),AL(170),BE(170),GA(170),P(170),Q(170),R(170),
     7  YC(170),SRESID(170),QCORCT(170)
C
      LOGICAL NVT,NVB
C
      EQUIVALENCE (RHS1(1),AJSQ(1)),(RHS2(1),AL(1)),(PCORCT(1),C(1)),
     1  (RHS3  (1),BE    (1),XF    (1)),(RHS4  (1),GA    (1),YG    (1)),
     2  (PRESID(1),F     (1),XXD   (1)),(QRESID(1),G     (1),YXD   (1)),
     3  (RRESID(1),A     (1),XED   (1)),(SRESID(1),B     (1),YED   (1)),
     4  (QCORCT(1),P     (1),XC    (1)),(RCORCT(1),Q     (1),YC    (1)),
     5  (SCORCT(1),R     (1),D     (1)),(S     (1),H     (1))
C
      DATA NCALL,NPANIC/0,0/
      DATA NLPP/60/
C
C        SET INDEX LIMITS, FACTORS FOR USE IN DIFFERENCE FORMULAS, ETC.
      NCALL=NCALL+1
      RDX=1./FLOAT(JINC)
      R2DX=.5*RDX
      RDXSQ=1./FLOAT(JINC)**2
      RDE=1./FLOAT(KINC)
      R2DE=.5*RDE
      RDESQ=1./FLOAT(KINC)**2
      RDXDE=1./(FLOAT(JINC)*FLOAT(KINC))
      TEORDX=TEOPEN*RDX
      TEO2DX=TEOPEN*R2DX
      TEODXS=TEOPEN*RDXSQ
      TALMR=TAN(ALAMR)
      CALMR=COS(ALAMR)
      SALMR=SIN(ALAMR)
      TERDXC=TEORDX*CALMR
      TERDXS=TEORDX*SALMR
      TE2DXS=TEO2DX*SALMR
      TE2DXC=TEO2DX*CALMR
      TEDXSC=TEODXS*CALMR
      TEDXSS=TEODXS*SALMR
      K2=KS+KINC
      K3=KS+2*KINC
      K4=KE-2*KINC
      K5=KE-KINC
      J2=JS+JINC
      J3=JE-JINC
      IF(NTETYP.NE.3.OR.NOBSHP.NE.6.AND.JTEBOT.NE.1)  GO TO 99
      NVB=(X(1,KMAX).NE.X(1,1))
      NVT=(X(JMAX,KMAX).NE.X(JMAX,1))
      IF(NVB)  SLABB=(Y(1,KMAX)-Y(1,1))/(X(1,KMAX)-X(1,1))
      IF(NVT)  SLABT=(Y(JMAX,KMAX)-Y(JMAX,1))/(X(JMAX,KMAX)-X(JMAX,1))
   99 CONTINUE
      JN1P=JN1+JINC
      JN1M=JN1-JINC
      JN2P=JN2+JINC
      ONE3RD=1./3.
      TWO3RD=2./3.
      NLPPM=NLPP-3
C
      DO 8 J=JS,JE,JINC
      JPPER(J)=J+JINC
    8 JMPER(J)=J-JINC
      IF(NTETYP.EQ.3)  GO TO 5
      JPPER(JE)=JS
      JMPER(JS)=JE
C
    5 IF(NTETYP.GT.1.OR.JINC.GT.1)  GO TO 74
      JPQF=J2
      JPQL=JE
      GO TO 12
   74 IF(NTETYP.EQ.3)  GO TO 75
      JPQF=J2
      JPQL=J3
      GO TO 12
   75 JPQF=JS
      JPQL=JE
C
C        COMPUTE DERIVATIVES OF X AND Y WITH RESPECT TO ETA AT INNER
C        AND OUTER BOUNDARIES.  THESE ARE FIXED FOR ALL ITERATIONS
C        AND HAVE EMBEDED IN THEM THE GIVEN REQUIREMENTS OF ANGLES AND
C        SPACING AT THOSE BOUNDARIES.
   12 DO 50 J=JS,JE,JINC
      JP=JPPER(J)
      JM=JMPER(J)
      XXB(J)=(X(JP,KS)-X(JM,KS))*R2DX
      YXB(J)=(Y(JP,KS)-Y(JM,KS))*R2DX
      XXC(J)=(X(JP,KE)-X(JM,KE))*R2DX
      YXC(J)=(Y(JP,KE)-Y(JM,KE))*R2DX
      XXXB(J)=(X(JP,KS)-2.*X(J,KS)+X(JM,KS))*RDXSQ
      YXXB(J)=(Y(JP,KS)-2.*Y(J,KS)+Y(JM,KS))*RDXSQ
      XXXC(J)=(X(JP,KE)-2.*X(J,KE)+X(JM,KE))*RDXSQ
      YXXC(J)=(Y(JP,KE)-2.*Y(J,KE)+Y(JM,KE))*RDXSQ
      IF(NTETYP.EQ.3.OR.J.NE.JS.AND.J.NE.JE)  GO TO 22
      XXB(J)=XXB(J)+TE2DXS
      YXB(J)=YXB(J)+TE2DXC
      XXC(J)=XXC(J)+TE2DXS
      YXC(J)=YXC(J)+TE2DXC
   20 IF(J.EQ.JE)  GO TO 21
      XXXB(J)=XXXB(J)-TEDXSS
      YXXB(J)=YXXB(J)-TEDXSC
      XXXC(J)=XXXC(J)-TEDXSS
      YXXC(J)=YXXC(J)-TEDXSC
      GO TO 22
   21 XXXB(J)=XXXB(J)+TEDXSS
      YXXB(J)=YXXB(J)+TEDXSC
      XXXC(J)=XXXC(J)+TEDXSS
      YXXC(J)=YXXC(J)+TEDXSC
   22 STH=SIN(THETA(J))
      CTH=COS(THETA(J))
      DUM=DS(J)/SQRT(XXB(J)**2+YXB(J)**2)
      XEB(J)=DUM*(-XXB(J)*CTH-YXB(J)*STH)
      YEB(J)=DUM*(-YXB(J)*CTH+XXB(J)*STH)
      STH=SIN(THETOB(J))
      CTH=COS(THETOB(J))
      DUM=DSOB(J)/SQRT(XXC(J)**2+YXC(J)**2)
      XEC(J)=DUM*(-XXC(J)*CTH-YXC(J)*STH)
      YEC(J)=DUM*(-YXC(J)*CTH+XXC(J)*STH)
   50 CONTINUE
C
C        INSERT REAR BOUNDARY CONDITIONS FOR THE CASE OF A C-GRID.
      IF(NTETYP.LT.3)  GO TO 13
      ATABB=ATAN2(X(1,KMAX)-X(1,1),Y(1,KMAX)-Y(1,1))
      SABB=SIN(ATABB)
      CABB=COS(ATABB)
      XEB(1)=-DS(1)*SABB
      YEB(1)=-DS(1)*CABB
      XEC(1)=DSOB(1)*SABB
      YEC(1)=DSOB(1)*CABB
      ATABT=ATAN2(X(JMAX,KMAX)-X(JMAX,1),Y(JMAX,KMAX)-Y(JMAX,1))
      SABT=SIN(ATABT)
      CABT=COS(ATABT)
      XEB(JMAX)=-DS(JMAX)*SABT
      YEB(JMAX)=-DS(JMAX)*CABT
      XEC(JMAX)=DSOB(JMAX)*SABT
      YEC(JMAX)=DSOB(JMAX)*CABT
C
C        DIFFERENCE WITH RESPECT TO XI TO GET MIXED SECOND PARTIALS
C        OF X AND Y AT INNER AND OUTER BOUNDARIES.  STORE FACTORS
C        FOR USE IN COMPUTING P1, Q1, R1, AND S1 FOR EACH ITERATION.
   13 DO 1 J=JS,JE,JINC
      JP=JPPER(J)
      JM=JMPER(J)
      XXEB=(XEB(JP)-XEB(JM))*R2DX
      YXEB=(YEB(JP)-YEB(JM))*R2DX
      XXEC=(XEC(JP)-XEC(JM))*R2DX
      YXEC=(YEC(JP)-YEC(JM))*R2DX
      RJB(J)=1.0/(XXB(J)*YEB(J)-XEB(J)*YXB(J))
      RJSQ=RJB(J)**2
      AD=XEB(J)**2+YEB(J)**2
      BD=XXB(J)*XEB(J)+YXB(J)*YEB(J)
      GD=XXB(J)**2+YXB(J)**2
      F0(J)=-GD*RJSQ
      F1(J)=(-AD*XXXB(J)+2.0*BD*XXEB)*RJSQ
      F2(J)=(-AD*YXXB(J)+2.0*BD*YXEB)*RJSQ
      RJC(J)=1.0/(XXC(J)*YEC(J)-XEC(J)*YXC(J))
      RJSQ=RJC(J)**2
      AD=XEC(J)**2+YEC(J)**2
      BD=XXC(J)*XEC(J)+YXC(J)*YEC(J)
      GD=XXC(J)**2+YXC(J)**2
      F3(J)=-GD*RJSQ
      F4(J)=(-AD*XXXC(J)+2.0*BD*XXEC)*RJSQ
      F5(J)=(-AD*YXXC(J)+2.0*BD*YXEC)*RJSQ
    1 CONTINUE
C
C        *************************************************************
C        *                                                           *
C        *              START OF THE MAIN ITERATION LOOP             *
C        *                                                           *
C        *************************************************************
C
C        INSIDE OF THE MAIN ITERATION LOOP THERE ARE MANY SMALL LOOPS
C        WITH J AS THE INDEX.  THE CODE WAS WRITTEN THIS WAY TO
C        ALLOW AS MANY AS POSSIBLE OF THOSE LOOPS TO BE INSTACK ON THE
C        CDC 7600.  ON OTHER COMPUTERS IT MIGHT BE DESIRABLE TO
C        SIMPLIFY THE CODE BY CONSOLIDATING MANY OF THOSE LOOPS.
   29 CONTINUE
      DO 23 ITER=1,MAXIT
      CMAX=0.
      CSUM=0.
      Q1SUM=0.0
      P1CM=0.
      Q1CM=0.
      R1CM=0.
      S1CM=0.
      P1M=0.
      Q1M=0.
      R1M=0.
      S1M=0.
C
C        COMPUTE SECOND DERIVATIVES OF X AND Y WITH RESPECT TO ETA,
C        FROM WHICH P1, Q1, R1, AND S1 AT THIS ITERATION ARE COMPUTED.
C        APPLY THE LIMITING ALGORITHM TO COMBAT THE TENDANCY TOWARD
C        INSTABILITY IN THE FIRST FEW ITERATIONS.  ALSO COMPUTE MAXIMA
C        OF THOSE FUNCTIONS FOR USE ON THE PRINTOUT.
      DO 63 J=JPQF,JPQL,JINC
      XEEB=.5*RDESQ*(-7.*X(J,KS)+8.*X(J,K2)-X(J,K3))-3.*RDE*XEB(J)
   63 RHS1(J)=F1(J)+F0(J)*XEEB
      DO 52 J=JPQF,JPQL,JINC
      YEEB=.5*RDESQ*(-7.*Y(J,KS)+8.*Y(J,K2)-Y(J,K3))-3.*RDE*YEB(J)
   52 RHS2(J)=F2(J)+F0(J)*YEEB
      DO 72 J=JPQF,JPQL,JINC
   72 PRESID(J)=(YEB(J)*RHS1(J)-XEB(J)*RHS2(J))*RJB(J)-P1(J)
      DO 30 J=JPQF,JPQL,JINC
   30 QRESID(J)=(-YXB(J)*RHS1(J)+XXB(J)*RHS2(J))*RJB(J)-Q1(J)
      DO 70 J=JPQF,JPQL,JINC
      PJUMP=PLIM*AMAX1(ABS(P1(J)),1.0)
   70 PCORCT(J)=SIGN(AMIN1(OMEGP*ABS(PRESID(J)),PJUMP),PRESID(J))
      DO 31 J=JPQF,JPQL,JINC
      QJUMP=QLIM*AMAX1(ABS(Q1(J)),1.0)
   31 QCORCT(J)=SIGN(AMIN1(OMEGQ*ABS(QRESID(J)),QJUMP),QRESID(J))
      DO 64 J=JS,JE,JINC
      XEEC=.5*RDESQ*(-7.*X(J,KE)+8.*X(J,K5)-X(J,K4))+3.*RDE*XEC(J)
   64 RHS3(J)=F4(J)+F3(J)*XEEC
      DO 32 J=JS,JE,JINC
      YEEC=.5*RDESQ*(-7.*Y(J,KE)+8.*Y(J,K5)-Y(J,K4))+3.*RDE*YEC(J)
   32 RHS4(J)=F5(J)+F3(J)*YEEC
      DO 73 J=JS,JE,JINC
   73 RRESID(J)=(YEC(J)*RHS3(J)-XEC(J)*RHS4(J))*RJC(J)-R1(J)
      DO 33 J=JS,JE,JINC
   33 SRESID(J)=(-YXC(J)*RHS3(J)+XXC(J)*RHS4(J))*RJC(J)-S1(J)
      DO 35 J=JS,JE,JINC
      RJUMP=RLIM*AMAX1(ABS(R1(J)),1.0)
   35 RCORCT(J)=SIGN(AMIN1(OMEGR*ABS(RRESID(J)),RJUMP),RRESID(J))
      DO 71 J=JS,JE,JINC
      SJUMP=SLIM*AMAX1(ABS(S1(J)),1.0)
   71 SCORCT(J)=SIGN(AMIN1(OMEGS*ABS(SRESID(J)),SJUMP),SRESID(J))
      DO 37 J=JPQF,JPQL,JINC
      P1(J)=P1(J)+PCORCT(J)
      IF(JPRT.LT.0)  GO TO 37
      ABSPC=ABS(PCORCT(J))
      IF(ABSPC.LT.P1CM)  GO TO 36
      P1CM=ABSPC
   36 ABSP1=ABS(P1(J))
      IF(ABSP1.LT.P1M)  GO TO 37
      P1M=ABSP1
      JP1M=J
   37 CONTINUE
      DO 39 J=JPQF,JPQL,JINC
      Q1(J)=Q1(J)+QCORCT(J)
      IF(JPRT.LT.0)  GO TO 39
      ABSQC=ABS(QCORCT(J))
      IF(ABSQC.LT.Q1CM)  GO TO 38
      Q1CM=ABSQC
   38 ABSQ1=ABS(Q1(J))
      IF(ABSQ1.LT.Q1M)  GO TO 39
      Q1M=ABSQ1
      JQ1M=J
   39 Q1SUM=Q1SUM+Q1(J)
      DO 41 J=JS,JE,JINC
      R1(J)=R1(J)+RCORCT(J)
      IF(JPRT.LT.0)  GO TO 41
      ABSRC=ABS(RCORCT(J))
      IF(ABSRC.LT.R1CM)  GO TO 40
      R1CM=ABSRC
   40 ABSR1=ABS(R1(J))
      IF(ABSR1.LT.R1M)  GO TO 41
      R1M=ABSR1
      JR1M=J
   41 CONTINUE
      DO 43 J=JS,JE,JINC
      S1(J)=S1(J)+SCORCT(J)
      IF(JPRT.LT.0)  GO TO 43
      ABSSC=ABS(SCORCT(J))
      IF(ABSSC.LT.S1CM)  GO TO 42
      S1CM=ABSSC
   42 ABSS1=ABS(S1(J))
      IF(ABSS1.LT.S1M)  GO TO 43
      S1M=ABSS1
      JS1M=J
   43 CONTINUE
C
C        IN THE CASE OF AN O-GRID WITH A SHARP TRAILING EDGE AND A
C        GRID POINT AT THAT TRAILING EDGE AND THE FINE-MESH SOLUTION,
C        COMPROMISE P1 AND Q1 AT THE TRAILING EDGE TO FACILITATE
C        CONVERGENCE.
      IF(NTETYP.NE.1.OR.JINC.NE.1)  GO TO 3
      P1(1)=0.5*(P1(2)+P1(JMAX))
      Q1(1)=0.5*(P1(2)+Q1(JMAX))
      GO TO 14
C
C        FOR ALL OTHER O-GRID CASES, COMPROMISE P1 AND Q1 AT THE
C        TRAILING EDGE TO FACILITATE CONVERGENCE.
    3 IF(NTETYP.EQ.3)  GO TO 14
      PDIF=P1(J2)-P1(J3)
      P1(JS)=P1(J3)+PDIF*TWO3RD
      P1(JE)=P1(J3)+PDIF*ONE3RD
      QDIF=Q1(J2)-Q1(J3)
      Q1(JS)=Q1(J3)+QDIF*TWO3RD
      Q1(JE)=Q1(J3)+QDIF*ONE3RD
C
C        COMPROMISE P1 AND Q1 AT LEADING EDGE TO FACILITATE CONVERGENCE
C        IF LEADING EDGE IS SHARP.
   14 GO TO (9,18,19),NLETYP
   18 P1(JN1)=.5*(P1(JN1M)+P1(JN1P))
      Q1(JN1)=.5*(Q1(JN1M)+Q1(JN1P))
      GO TO 9
   19 PDIF=P1(JN2P)-P1(JN1M)
      P1(JN1)=P1(JN1M)+PDIF*ONE3RD
      P1(JN2)=P1(JN1M)+PDIF*TWO3RD
      QDIF=Q1(JN2P)-Q1(JN1M)
      Q1(JN1)=Q1(JN1M)+QDIF*ONE3RD
      Q1(JN2)=Q1(JN1M)+QDIF*TWO3RD
    9 CONTINUE
C
C        LOOP ON K.  DIRECTION VARRIES WITH SIGN OF Q1SUM.
      DO 10 KK=K2,K5,KINC
      IF(Q1SUM.GT.0.0)  GO TO 26
      K=KK
      GO TO 27
   26 K=KE+KS-KK
   27 KP=K+KINC
      KM=K-KINC
C
C        FILL TRIDIAGONAL ARRAYS.
      DO 65 J=JS,JE,JINC
      JP=JPPER(J)
      JM=JMPER(J)
      XXD(J)=(X(JP,K)-X(JM,K))*R2DX
   65 YXD(J)=(Y(JP,K)-Y(JM,K))*R2DX
      XXD(JS)=XXD(JS)+TE2DXS
      YXD(JS)=YXD(JS)+TE2DXC
      XXD(JE)=XXD(JE)+TE2DXS
      YXD(JE)=YXD(JE)+TE2DXC
      DO 11 J=JS,JE,JINC
      XED(J)=(X(J,KP)-X(J,KM))*R2DE
   11 YED(J)=(Y(J,KP)-Y(J,KM))*R2DE
      DO 45 J=JS,JE,JINC
      AJSQ(J)=(XXD(J)*YED(J)-XED(J)*YXD(J))**2
      AL(J)=XED(J)**2+YED(J)**2
      BE(J)=XXD(J)*XED(J)+YXD(J)*YED(J)
   45 GA(J)=XXD(J)**2+YXD(J)**2
      DO 66 J=JS,JE,JINC
      JP=JPPER(J)
      JM=JMPER(J)
      XXED=.25*RDXDE*(X(JP,KP)-X(JP,KM)-X(JM,KP)+X(JM,KM))
   66 F(J)=2.*BE(J)*XXED
      DO 46 J=JS,JE,JINC
      JP=JPPER(J)
      JM=JMPER(J)
      YXED=.25*RDXDE*(Y(JP,KP)-Y(JP,KM)-Y(JM,KP)+Y(JM,KM))
   46 G(J)=2.*BE(J)*YXED
      DO 69 J=JS,JE,JINC
      F(J)=F(J)-GA(J)*RDESQ*(X(J,KP)+X(J,KM))
   69 G(J)=G(J)-GA(J)*RDESQ*(Y(J,KP)+Y(J,KM))
      DO 47 J=JS,JE,JINC
      A(J)=AL(J)*RDXSQ
      B(J)=-2.*AL(J)*RDXSQ-2.*GA(J)*RDESQ
   47 C(J)=AL(J)*RDXSQ
      DO 48 J=JS,JE,JINC
      P(J)=P1(J)*EXP(-(K-1)*AAA(J))
      Q(J)=Q1(J)*EXP(-(K-1)*BBB(J))
      R(J)=R1(J)*EXP((K-KE)*CCC(J))
   48 S(J)=S1(J)*EXP((K-KE)*DDD(J))
      DO 49 J=JS,JE,JINC
      PPR=P(J)+R(J)
      ADD=AJSQ(J)*ABS(PPR)*RDX
      B(J)=B(J)-ADD
      IF(PPR.GE.0.)  GO TO 62
      A(J)=A(J)+ADD
      GO TO 49
   62 C(J)=C(J)+ADD
   49 CONTINUE
      DO 25 J=JS,JE,JINC
      QPS(J)=Q(J)+S(J)
      IF(QPS(J).GT.0.0)  GO TO 24
      XF(J)=-X(J,KM)
      GO TO 25
   24 XF(J)=X(J,KP)
   25 CONTINUE
      DO 67 J=JS,JE,JINC
      IF(QPS(J).GT.0.0)  GO TO 68
      YG(J)=-Y(J,KM)
      GO TO 67
   68 YG(J)=Y(J,KP)
   67 CONTINUE
      DO 60 J=JS,JE,JINC
      B(J)=B(J)-AJSQ(J)*ABS(QPS(J))*RDE
      F(J)=F(J)-AJSQ(J)*QPS(J)*XF(J)*RDE
   60 G(J)=G(J)-AJSQ(J)*QPS(J)*YG(J)*RDE
C
C        IN THE CASE OF C-GRIDS, INSERT REAR BOUNDARY CONDITIONS
C        AND SOLVE TRIDIAGONAL SYSTEMS.
      IF(NTETYP.LT.3)  GO TO 15
      F(JS)=F(JS)-A(JS)*X(1,K)
      F(JE)=F(JE)-C(JE)*X(JMAX,K)
      G(1)=(X(1,K)-X(JS,K))*TALMR
      G(JMAX)=(X(JMAX,K)-X(JE,K))*TALMR
      A(1)=0.
      B(1)=-1.
      C(1)=1.
      A(JMAX)=1.
      B(JMAX)=-1.
      C(JMAX)=0.
      CALL TRIB(A,B,C,D,F,JS,JE,JINC)
      CALL TRIB(A,B,C,D,G,1,JMAX,JINC)
      GO TO 16
C
C        SOLVE TRIDIAGONAL SYSTEMS WITH PERIODIC BOUNDARY CONDITIONS
C        FOR O-GRIDS.
   15 CONTINUE
      F(JS)=F(JS)+AL(JS)*TEDXSS
      G(JS)=G(JS)+AL(JS)*TEDXSC
      F(JE)=F(JE)-AL(JE)*TEDXSS
      G(JE)=G(JE)-AL(JE)*TEDXSC
      PPR=P(JS)+R(JS)
      IF(PPR.LT.0.)  F(JS)=F(JS)-AJSQ(JS)*PPR*TERDXS
      IF(PPR.LT.0.)  G(JS)=G(JS)-AJSQ(JS)*PPR*TERDXC
      PPR=P(JE)+R(JE)
      IF(PPR.GE.0.)  F(JE)=F(JE)-AJSQ(JE)*PPR*TERDXS
      IF(PPR.GE.0.)  G(JE)=G(JE)-AJSQ(JE)*PPR*TERDXC
      CALL TRIP(A,B,C,F,D,H,JS,JE,JINC)
      CALL TRIP(A,B,C,G,D,H,JS,JE,JINC)
C
C        UPDATE X AND Y ARRAYS AT THIS K
   16 DO 61 J=JS,JE,JINC
      XC(J)=F(J)-X(J,K)
      YC(J)=G(J)-Y(J,K)
      X(J,K) = X(J,K) + OMEGA*XC(J)
   61 Y(J,K) = Y(J,K) + OMEGA*YC(J)
C
C        FOR C-TYPE GRIDS INSERT REAR BOUNDARY CONDITIONS
      IF(NTETYP.LT.3)  GO TO 17
      IF(NOBSHP.EQ.6.OR.JTEBOT.EQ.1)  GO TO 34
      Y(1,K)=G(1)
      Y(JMAX,K)=G(JMAX)
      GO TO 17
   34 SL12B=(Y(JS,K)-Y(J2,K))/(X(JS,K)-X(J2,K))
      SL12T=(Y(JE,K)-Y(J3,K))/(X(JE,K)-X(J3,K))
      IF(NVB)  X(1,K)=(X(J2,K)*SL12B-X(1,1)*SLABB
     1  +Y(1,1)-Y(J2,K))/(SL12B-SLABB)
      IF(NVT)  X(JMAX,K)=(X(J3,K)*SL12T
     1  -X(JMAX,1)*SLABT+Y(JMAX,1)-Y(J3,K))/(SL12T-SLABT)
      Y(1,K)=SL12B*(X(1,K)-X(J2,K))+Y(J2,K)
      Y(JMAX,K)=SL12T*(X(JMAX,K)-X(J3,K))+Y(J3,K)
   17 CONTINUE
C
C        COMPUTE SUM OF CORRECTIONS.
      DO 4 J=JS,JE,JINC
      CORECT=ABS(XC(J))+ABS(YC(J))
      IF(CORECT.LE.CMAX)  GO TO 4
      CMAX=CORECT
      JCMAX=J
      KCMAX=K
    4 CSUM=CSUM+CORECT
   10 CONTINUE
C        END OF LOOP ON K.
C
C        PRINT FOR THIS ITERATION
      IF(JPRT.LT.0)  GO TO 6
      IF((ITER-1)/NLPPM*NLPPM.NE.ITER-1)  GO TO 44
      IF(JINC.NE.1.OR.KINC.NE.1)  WRITE(6,111)
      IF(JINC.EQ.1.AND.KINC.EQ.1)  WRITE(6,104)
   44 CONTINUE
      WRITE(6,100)  ITER,CSUM,CMAX,JCMAX,KCMAX,P1CM,P1M,JP1M,
     1  Q1CM,Q1M,JQ1M,R1CM,R1M,JR1M,S1CM,S1M,JS1M
    6 CONTINUE
C
C        IF MAXIMUM CORRECTION ON X AND Y IS REDUCED NORDA(NCALL)
C        ORDERS OF MAGNITUDE, QUIT.
      IF(ITER.GT.1)  GO TO 28
      CMSTRT=CMAX
      CMQUIT=CMAX*10.**(-NORD)
      CPANIC=CMAX*10.
      GO TO 23
   28 IF(CMAX.LE.CMQUIT)  GO TO 2
C
C        RESTART THE SOLUTION IF IT IS GOING UNSTABLE.
      IF(CMAX.LT.CPANIC.OR.NCALL.EQ.2.OR.NPANIC.EQ.3)  GO TO 23
      IF(JPRT.GE.0)  WRITE(6,108)
      NPANIC=NPANIC+1
      PLIM=PLIM*.5
      QLIM=QLIM*.5
      RLIM=RLIM*.5
      SLIM=SLIM*.5
      OMEGA=OMEGA*.75
      OMEGP=OMEGP*.75
      OMEGQ=OMEGQ*.75
      OMEGR=OMEGR*.75
      OMEGS=OMEGS*.75
      IF(JPRT.GE.0)  WRITE(6,109)  NPANIC,PLIM,QLIM,RLIM,SLIM,OMEGA,
     1  OMEGP,OMEGQ,OMEGR,OMEGS
      CALL IC
      GO TO 29
C
   23 CONTINUE
      ITER=MAXIT+1
C
C        *************************************************************
C        *                                                           *
C        *                END OF THE MAIN ITERATION LOOP             *
C        *                                                           *
C        *************************************************************
C
C        PRINT ONE-LINE CONVERGENCE MESSAGE.
    2 CONTINUE
      IF(ITER.LE.MAXIT)  GO TO 80
      IF(JINC.EQ.1.AND.KINC.EQ.1)  GO TO 82
      NCMRED=IFIX(ALOG10(CMSTRT/CMAX)+.5)
      NDIF=(NORD-NCMRED)/2
      NORDA(2)=NORDA(2)+NDIF
      WRITE(6,110)  MAXIT,NCMRED,NORD,NDIF,NORDA(2)
      GO TO 81
   82 WRITE(6,103)  NORD,MAXIT
      GO TO 81
   80 IF(JINC.NE.1.OR.KINC.NE.1)  WRITE(6,112)  NORD,ITER
      IF(JINC.EQ.1.AND.KINC.EQ.1)  WRITE(6,107)  NORD,ITER
   81 CONTINUE
C
C        PRINT ANGLES, ETC., OF SOLUTION.
      IF(JPRT.LT.0)  GO TO 7
      N=0
C
      DO 97 J=JS,JE,JINC
      N=N+1
      IF((N-1)/NLPPM*NLPPM.EQ.N-1)  WRITE(6,105)  KS,KS,K2,K2
      JP=JPPER(J)
      JM=JMPER(J)
      ANG1=ATAN2(Y(J,K2)-Y(J,KS),X(J,K2)-X(J,KS))*TODEG
      IF(ANG1.LT.0.0)  ANG1=ANG1+360.0
      ANG2=ATAN2(Y(JP,KS)-Y(JM,KS),X(JP,KS)-X(JM,KS))*TODEG
      IF(ANG2.LT.0.0)  ANG2=ANG2+360.0
      ANGA=ANG2-ANG1+180.
      IF(ANGA.GT.360.)  ANGA=ANGA-360.
      IF(ANGA.LT.0.0)  ANGA=ANGA+360.
      DELS=SQRT((X(J,K2)-X(J,KS))**2+(Y(J,K2)-Y(J,KS))**2)
   97 WRITE(6,101)  J,ANGA,ANG1,DELS,X(J,KS),Y(J,KS),X(J,K2),Y(J,K2),
     1  P1(J),Q1(J)
      N=0
C
      DO 98 J=JS,JE,JINC
      N=N+1
      IF((N-1)/NLPPM*NLPPM.EQ.N-1)  WRITE(6,106)  KE,KE,K5,K5
      JP=JPPER(J)
      JM=JMPER(J)
      ANG1=ATAN2(Y(J,KE)-Y(J,K5),X(J,KE)-X(J,K5))*TODEG
      IF(ANG1.LT.0.0)  ANG1=ANG1+360.
      ANG2=ATAN2(Y(JP,KE)-Y(JM,KE),X(JP,KE)-X(JM,KE))*TODEG
      IF(ANG2.LT.0.0)  ANG2=ANG2+360.
      ANGA=ANG1-ANG2
      IF(ANGA.LT.0.0)  ANGA=ANGA+360.
      ANGA=180.-ANGA
      DELS=SQRT((X(J,K5)-X(J,KE))**2+(Y(J,K5)-Y(J,KE))**2)
   98 WRITE(6,101)  J,ANGA,ANG1,DELS,X(J,KE),Y(J,KE),X(J,K5),Y(J,K5),
     1  R1(J),S1(J)
    7 CONTINUE
C
      RETURN
C
  100 FORMAT(1X,I4,1P2E10.3,I4,I3,4(1X,2E10.3,I4))
  101 FORMAT(1X,I3,2F12.2,F12.8,4F12.6,2E15.6)
  103 FORMAT(29H0FINE SOLUTION NOT CONVERGED./24H MAXIMUM CORRECTION NOT
     * ,
     1  21HREDUCED THE REQUIRED ,I2,24H ORDERS OF MAGNITUDE IN ,I5,
     2  12H ITERATIONS.)
  104 FORMAT(1H1,T40,26HFINE GRID SOLUTION PROCESS/
     3  5H ITER,6X,4HCSUM,6X,4HCMAX,3X,1HJ,2X,1HK,7X,4HP1CM,
     1  7X,3HP1M,3X,1HJ,7X,4HQ1CM,7X,3HQ1M,3X,1HJ,7X,4HR1CM,7X,3HR1M,
     2  3X,1HJ,7X,4HS1CM,7X,3HS1M,3X,1HJ)
  105 FORMAT(1H1,50X,37HSOLUTION AT INNER (AIRFOIL) BOUNDARY:/
     1  3X,1HJ,3X,9HANGLE W/R,3X,9HANGLE W/R,12H DISTANCE TO,
     2  5X,4HX(J,,I2,1H),5X,4HY(J,,I2,1H),5X,4HX(J,,I2,1H),5X,4HY(J,,
     3  I2,1H),13X,2HP1,13X,2HQ1/8X,8HBOUNDARY,6X,6HX-AXIS,4X,
     4  8HBOUNDARY)
  106 FORMAT(1H1,50X,27HSOLUTION AT OUTER BOUNDARY:/
     1  3X,1HJ,3X,9HANGLE W/R,3X,9HANGLE W/R,12H DISTANCE TO,
     2  5X,4HX(J,,I2,1H),5X,4HY(J,,I2,1H),5X,4HX(J,,I2,1H),5X,4HY(J,,
     3  I2,1H),13X,2HR1,13X,2HS1/8X,8HBOUNDARY,6X,6HX-AXIS,4X,
     4  8HBOUNDARY)
  107 FORMAT(54H0FINE SOLUTION CONVERGED.  MAXIMUM CORRECTION REDUCED ,
     1 I2,24H ORDERS OF MAGNITUDE IN ,I5,12H ITERATIONS.)
  108 FORMAT(////47H SOLUTION PROCESS APPEARS TO BE GOING UNSTABLE./
     1  39H RELAXATION PARAMETERS WILL BE REDUCED./
     2  42H GRID WILL BE RESET TO INITIAL CONDITOINS./
     3  37H ITERATION PROCESS WILL BE RESTARTED.)
  109 FORMAT(8H0NPANIC=,I2,5X,5HPLIM=,F7.3,5X,5HQLIM=,F7.3,5X,5HRLIM=,
     1  F7.3,5X,5HSLIM=,F7.3/7H0OMEGA=,F7.3,5X,6HOMEGP=,F7.3,5X,
     2  6HOMEGQ=,F7.3,5X,6HOMEGR=,F7.3,5X,6HOMEGS=,F7.3)
  110 FORMAT(34H0COARSE SOLUTION NOT CONVERGED IN ,I5,14H ITERATIONS.  ,
     Z  7HMAXIMUM,
     1  27H CORRECTION REDUCED APROX. ,I2,20H ORDERS OF MAGNITUDE/
     2  25H INSTEAD OF THE REQUIRED ,I2,27H.  HALF OF THE DIFFERENCE, ,I
     *2,
     3  46H, HAS BEEN ADDED TO NORDA(2) WHICH NOW EQUALS ,I2,1H.)
  111 FORMAT(1H1,T40,28HCOARSE GRID SOLUTION PROCESS/
     3  5H ITER,6X,4HCSUM,6X,4HCMAX,3X,1HJ,2X,1HK,7X,4HP1CM,
     1  7X,3HP1M,3X,1HJ,7X,4HQ1CM,7X,3HQ1M,3X,1HJ,7X,4HR1CM,7X,3HR1M,
     2  3X,1HJ,7X,4HS1CM,7X,3HS1M,3X,1HJ)
  112 FORMAT(56H0COARSE SOLUTION CONVERGED.  MAXIMUM CORRECTION REDUCED
     *,
     1 I2,24H ORDERS OF MAGNITUDE IN ,I5,12H ITERATIONS.)
      END
*DECK SOLVE
      SUBROUTINE SOLVE
C
C        IF GRID SEQUENCING IS TO BE USED, THIS SUBROUTINE SUPPLIES
C        INITIAL CONDITIONS FOR X AND Y BY CALLING IC, THEN SOLVES
C        THE EQUATIONS ON THE COARSE GRID BY CALLING RELAX FOR THE
C        FIRST TIME, THEN INTERPOLATES THE SOLUTION TO FIND INITIAL
C        CONDITIONS FOR THE FINE GRID SOLUTION, THEN CALLS RELAX
C        AGAIN TO EFFECT THAT FINE GRID SOLUTION.  IF GRID SEQUENCING IS
C        NOT TO BE USED, THIS SUBROUTINE CALLS IC, THEN RELAX FOR THE
C        FINE GRID SOLUTION.
C
*CALL GRIDB
C
*CALL GRIDC
C
      KS=1
      KE=KMAX
C
C        SOLVE FOR THE GRID IN THE CASE OF AN O-GRID WITH A SHARP
C        TRAILING EDGE, AND A POINT AT THAT TRAILING EDGE.
      IF(NTETYP.GT.1)  GO TO 1
      IF(NSEQ.EQ.1)  GO TO 5
      JS=3
      JE=JMAX-1
      JINC=4
      KINC=3
      JN2=JLE2A(1)
      JN1=JLE1A(1)
      MAXIT=MAXITA(1)
      NORD=NORDA(1)
      CALL IC
      CALL RELAX
      CALL INTERP
    5 JS=1
      JE=JMAX
      JINC=1
      KINC=1
      JN1=JLE1A(2)
      JN2=JLE2A(2)
      MAXIT=MAXITA(2)
      NORD=NORDA(2)
      IF(NSEQ.EQ.1)  CALL IC
      IF(MAXIT.EQ.0)  GO TO 4
      CALL RELAX
      GO TO 4
C
C        SOLVE FOR THE GRID IN ALL OTHER O-GRID CASES.
    1 IF(NTETYP.EQ.3)  GO TO 7
      IF(NSEQ.EQ.1)  GO TO 8
      JS=2
      JE=JMAX-1
      JINC=3
      KINC=3
      JN2=JLE2A(1)
      JN1=JLE1A(1)
      MAXIT=MAXITA(1)
      NORD=NORDA(1)
      CALL IC
      CALL RELAX
      CALL INTERP
    8 JS=1
      JE=JMAX
      JINC=1
      KINC=1
      JN1=JLE1A(2)
      JN2=JLE2A(2)
      MAXIT=MAXITA(2)
      NORD=NORDA(2)
      IF(NSEQ.EQ.1)  CALL IC
      IF(MAXIT.EQ.0)  GO TO 4
      CALL RELAX
      GO TO 4
C
C        SOLVE FOR THE GRID IN THE CASE OF A C-GRID.
    7 IF(NSEQ.EQ.1)  GO TO 10
      JS=4
      JE=JMAX-3
      JINC=3
      KINC=3
      JN2=JLE2A(1)
      JN1=JLE1A(1)
      MAXIT=MAXITA(1)
      NORD=NORDA(1)
      CALL IC
      CALL RELAX
      CALL INTERP
   10 JS=2
      JE=JMAX-1
      JINC=1
      KINC=1
      JN1=JLE1A(2)
      JN2=JLE2A(2)
      MAXIT=MAXITA(2)
      NORD=NORDA(2)
      IF(NSEQ.EQ.1)  CALL IC
      IF(MAXIT.EQ.0)  GO TO 4
      CALL RELAX
    4 RETURN
C
      END
*DECK TRIB
      SUBROUTINE TRIB(A,B,C,X,F,NL,NU,JINC)
      DIMENSION A(2),B(2),C(2),X(2),F(2)
C
C        THIS SUBROUTINE SOLVES A TRI-DIAGONAL SYSTEM OF LINEAR
C        EQUATIONS.
C
      X(NL)=C(NL)/B(NL)
      F(NL)=F(NL)/B(NL)
      NLP1=NL+JINC
      DO 1 J=NLP1,NU,JINC
      JM=J-JINC
      Z=1./(B(J)-A(J)*X(JM))
      X(J)=C(J)*Z
    1 F(J)=(F(J)-A(J)*F(JM))*Z
      NUPNL=NU+NL
      DO 2 J1=NLP1,NU,JINC
      J=NUPNL-J1
      JP=J+JINC
    2 F(J)=F(J)-X(J)*F(JP)
      RETURN
      END
*DECK TRIP
      SUBROUTINE TRIP(A,B,C,F,Q,S,J1,J2,JINC)
      DIMENSION A(3),B(3),C(3),F(3),Q(3),S(3)
      JA=J1+JINC
      FN = F(J2)
C   FORWARD ELIMINATION SWEEP
      Q(J1) = -C(J1)/B(J1)
      F(J1) = F(J1)/B(J1)
      S(J1) = - A(J1)/B(J1)
      DO 10 J=JA,J2,JINC
      JM=J-JINC
      P=1./(B(J)+A(J)*Q(JM))
      Q(J) = - C(J)*P
      F(J)=(F(J)-A(J)*F(JM))*P
      S(J)=-A(J)*S(JM)*P
   10 CONTINUE
C   BACKWARD PASS
      JJ = J1 + J2
      Q(J2) = 0.
      S(J2) = 1.
      DO 11 I=JA,J2,JINC
      J = JJ - I
      JP=J+JINC
      S(J)=S(J)+Q(J)*S(JP)
   11 Q(J)=F(J)+Q(J)*Q(JP)
      J2M=J2-JINC
      F(J2)=(FN-C(J2)*Q(J1)-A(J2)*Q(J2M))/(C(J2)*S(J1)+A(J2)*S(J2M)
     1  +B(J2))
C    BACKWARD ELIMINATION PASS
      DO 12 I=JA,J2,JINC
      J = JJ -I
   12 F(J) = F(J2)*S(J) + Q(J)
      RETURN
      END
*DECK MAIN20
      OVERLAY(2,0)
C
      PROGRAM MAIN20
C
*CALL GRIDA
C
*CALL GRIDB
C
*CALL GRIDC
C
      IF(JPRT.GE.0)  WRITE(6,108)
      DO 8 N=1,NPLT
      CALL PLAWT(JMAX,KMAX,X,Y,XMIN(N),XMAX(N),YMIN(N),YMAX(N))
      IF(JPRT.GE.0)  WRITE(6,109)  N,XMIN(N),XMAX(N),YMIN(N),YMAX(N)
    8 CONTINUE
      CALL DONEPL
      RETURN
C
  108 FORMAT(14H0LIST OF PLOTS/10H  PLOT NO.,8X,4HXMIN,8X,4HXMAX,
     1  18X,4HYMIN,8X,4HYMAX)
  109 FORMAT(10X,I3,2F12.6,10X,2F12.6)
      END
*DECK PLAWT
      SUBROUTINE PLAWT(JMAX,KMAX,X,Y,XMIN,XMAX,YMIN,YMAX)
C
C        PLOT THE GRID USING ISSCO DISSPLA SOFTWARE.
C
      DIMENSION X(170,85),Y(170,85)
C
      DIMENSION XX(85),YY(85)
C
      DATA NPLT/0/
C
C        READJUST PLOT LIMITS TO AVIOD A STRECHED PLOT.
      XDIF=XMAX-XMIN
      YDIF=YMAX-YMIN
      IF(XDIF.LT.YDIF)  GO TO 4
      XDIFH=XDIF*0.5
      YMID=(YMAX+YMIN)*0.5
      YMX=YMID+XDIFH
      YMN=YMID-XDIFH
      XMX=XMAX
      XMN=XMIN
      GO TO 5
    4 YDIFH=YDIF*0.5
      XMID=(XMAX+XMIN)*0.5
      XMX=XMID+YDIFH
      XMN=XMID-YDIFH
      YMX=YMAX
      YMN=YMIN
    5 CONTINUE
C
C        PLOT THE LINES.
      NPLT=NPLT+1
      IF(NPLT.GT.1)  GO TO 6
      CALL BETA
      CALL BGNPL(-1)
      CALL NOCHEK
      CALL GRACE(0.001)
    6 CALL TITLE(1H ,1,0,0,0,0,7.5,7.5)
      CALL FRAME
      XYSCAL=(XMX-XMN)/7.5
      CALL GRAPH(XMN,XYSCAL,YMN,XYSCAL)
      DO 1 K=1,KMAX
    1 CALL CURVE(X(1,K),Y(1,K),JMAX,0)
      DO 2 J=1,JMAX
      DO 3 K=1,KMAX
      XX(K)=X(J,K)
    3 YY(K)=Y(J,K)
    2 CALL CURVE(XX,YY,KMAX,0)
C
C        ANNOTATE THE PLOT
      CALL MESSAG(4HXMN=,4,0.0,-0.75)
      CALL REALNO(XMN,2,4HABUT,4HABUT)
      CALL MESSAG(4HXMX=,4,2.0,-0.75)
      CALL REALNO(XMX,2,4HABUT,4HABUT)
      CALL MESSAG(4HYMN=,4,4.0,-0.75)
      CALL REALNO(YMN,2,4HABUT,4HABUT)
      CALL MESSAG(4HYMX=,4,6.0,-0.75)
      CALL REALNO(YMX,2,4HABUT,4HABUT)
      CALL ENDPL(-NPLT)
      RETURN
      END
